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ABSTRACT 

Integral  equations  (lE's)  are  widely  utilized  to  calculate  induced  currents  on  anten- 
nas and  scatterers,  but  they  are  seriously  restricted  in  their  ability  to  handle  inhomoge- 
neous  penetrable  structures  having  multiwavelength  dimensions.  The  utilization  of  finite 
element  (FE)  techniques  has  not  been  as  pervasive  as  the  use  of  lE's.  The  IE  represen- 
tation matrix  is  'Tull",  containing  few,  if  any,  zero  valued  elements.  The  techniques  for 
operating  on  these  large-sized  full  matrices  require  undesirable  amounts  of  processor 
time.  FE  techniques  produce  sparse  matrices  due  to  the  strictly  local  interactions  be- 
tween discrete  unknowns.  The  application  of  FE's  to  unbounded  problems,  however, 
requires  supplementary  enforcement  of  the  far-field  radiation  conditions.  The  Field 
Feedback  Formulation  {F  )  circumvents  the  full-matrix  computational  "bottleneck"  by 
allowing  FE  based  numerical  methods  to  be  employed.  Even  though  the  resultant  sparse 
matrices  may  be  larger  than  the  "full"  matrices  discussed  earlier,  most  elements  have  a 
value  of  zero.  Numerical  procedures  exist  to  optimize  operations  with  these  sparse  ma- 
trices. Calculational  speeds  can  be  orders  of  magnitude  faster.  Computer  techniques  to 
implement  and  validate  this  new  technique  are  the  basis  for  this  thesis.  Excellent 
agreement  with  classical  results  are  demonstrated. 
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I.     INTRODUCTION 

A.  HISTORY 

The  application  of  finite  element  techniques  to  evaluate  the  solution  of  differential 
equations  is  well  documented.  The  utilization  of  these  techniques  by  the 
electromagnetics  community  has  not  been  as  pervasive  as  the  use  of  integral  equations 
(lE's).  lE's  are  widely  utilized  to  calculate  induced  currents  on  antennas  and  scatterers, 
but  they  are  seriously  restricted  in  their  ability  to  handle  inhomogeneous  penetrable 
structures  having  multiwavelength  dimensions.  As  the  complexity  and  number  of  nodal 
degrees  of  freedom  grow,  the  size  and  dimension  of  the  representative  matrix  must  also 
grow.  This  matrix  is  "full",  containing  few,  if  any.  zero-valued  elements.  The  available 
numerical  techniques  for  operating  on  these  large-sized  full  matrices  require  undesirable 
amounts  of  processor  time. 

Differential  equation  (DE)  based  techniques,  such  as  the  finite  element  method, 
produce  sparse  matrices  due  to  the  strictly  local  interactions  between  discrete  unknowns 
which  result.  The  application  of  DE's  to  unbounded  problems,  such  as  those  of  scat- 
tering and  radiation,  require  some  form  of  supplementary  enforcement  of  the  proper 
far-field  conditions.  These  radiation  boundary  conditions  are  innately  incorporated  into 
integral  equations. 

B.  FIELD  FEEDBACK  FORMULATION 

The  Field  Feedback  Formulation  {F  )  circumvents  the  full-matrix  "bottleneck"  in  the 
computational  process  by  allowing  DE  based  numerical  methods  to  be  employed.  Even 
though  the  resultant  sparse  matrices  may  be  larger  than  the  "full"  matrices  discussed 
earlier,  most  elements  have  a  value  of  zero.  Numerical  procedures  exist  to  optimize 
operations  with  these  sparse  matrices.  Calculational  speeds  can  be  orders  of  magnitude 
faster  than  with  full  matrices  [Ref  1].  Although  the  F  employs  sparse  matrices  to  rep- 
resent the  fields  in  the  materials  being  considered,  it  does  require  augmentation  to  en- 
force the  radiation  condition  at  infinity  on  the  scattered  fields.  This  comes  in  the  form 
of  a  feedback  matrix  composed  of  surface  integration  generated  elements.  A  concept 
evaluation,  for  a  special  axisymmetric  case,  was  already  accomplished,  as  detailed  in 
[Ref  2]  and  [Ref  3].  Computer  techniques  to  implement  and  validate  this  new  technique 
are  the  basis  for  this  thesis. 


C.     POTENTIAL  BENEFITS 

This  thesis  will  lead  to  an  increased  understanding  of  the  advantages  and  disadvan- 
tages of  this  novel  computational  procedure  for  handling  geometrically  complex  material 
scatterers.  This  method  may  ultimately  allow  computer-aided  design  of  important 
electromagnetic  structures  such  as  low-observable  aircraft,  high  efficiency  dielectric  lens 
antennas,  and  other  electromagnetic  scattering  occurrences  due  to  atmospheric  anoma- 
lies. Structural  details  and  material  inhomogeneities,  as  well  as  physical  dimensions  (in 
muhiple  wavelengths),  can  be  accommodated  using  the  Field  Feedback  Formulation. 
These  capabilities  far  surpass  those  which  are  possible  with  contemporary  integral 
equation  techniques  for  the  case  of  inhomogeneous  penetrable  scatterers  and  antennas. 


II.     FORMULATION 

A.     INITIAL  NOMENCLATURE 

Assume  there  is  a  three  dimensional  object  that  is  infinite  in  one  direction.  Such  an 
object,  in  cross  section  could  look.  like  Figure  1.  This  object  only  varies  in  two  dimen- 
sions, and  therefore,  is  actually  a  two  dimensional  (2-D)  object. 
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Figure  1.      A  Typical  Object 


The  wavenumber  L  is  defined  as 


k  -^ 


27i/o 


where  }.„  is  the  fi-ee  space  wavelength  associated  with  an  electromagnetic  wave  of  fi^e- 
quency  f^  and  c  is  the  speed  of  light.     The  X  and  Y  Cartesian  coordinates  are 


wavenumber  normalized,  such  that  X  =  k^  and  Y  =  k^K  This  coordinate  normalization 
will  be  used  throughout  this  development.  A  similar  technique  could  also  be  developed 
in  the  polar  coordinate  system.  The  magnetic  field  will  also  be  normalized  such  that 
//  =  —j>]q-^,  where  >/(,  =  yj ~f'  =  \2Qn  x  311Q,  is  the  impedance  of  free  space  and 
y/f  is  the  usual  magnetic  field  in  units  of  A/m.  Thus  the  normalized  //  has  the  same  V/m 
units  as  E.   Potentials  may  be  defined  as, 


for  the  transverse  magnetic  (TM)  case,  with  £"„  E^,  and  //,  =  0  ,  and 


(2.1) 


(2.2) 


for  the  transverse  electric  (TE)  case,  with  II„  H^  and  E,  =  0  . 

Our  objective  is  to  calculate  the  scattered  fields  for  an  arbitrary  (2-D)  penetrable 
object  using  the  Field  Feedback  Formulation  (F  ).  As  shown  in  Figure  2,  a  familiar 
closed  loop  system  illustrates  the  relationship  between  the  incident,  scattered,  total  and 
far  fields. 
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Figure  2.      Field  Feedback  Formulation 


At  point  1,  the  incident  field  drives  the  total  system.  The  incident  field  at  1,  combined 
with  the  scattered  field  at  4,  forms  the  total  field  at  2.  This  total  field  forms  the 
boundary  conditions  that  drive  the  finite  element  program  at  2.  At  point  2,  initially,  the 

3 

incident  field  drives  the  U  operator.  U  represents  the  feed  forward  operator  in  the  F . 
This  operator  uses  a  finite  element  technique  to  solve  the  boundary  value  problem.  At 
point  3,  the  boundary  conditions  are  solved  for  the  object  perimeter  potentials  and  the 


normal  derivative  of  these  potentials.  T  represents  the  field  feedback  operator  that  takes 
the  perimeter  potentials  and  associated  derivatives  and  provides  the  scattered  fields  at 
point  4,  the  ofi^set  boundary.  The  fields  at  point  I  and  4  are  added  (unlike  the  familiar 
feedback  or  control  system  where  negative  feedback  is  employed).  At  point  2,  a  com- 
bined incident  and  scattered  field  exists.  These  fields  are  the  combined  boundary  con- 
ditions for  the  finite  element  boundan.'  value  program.  These  combined  fields  (total 
fields)  are  then  applied  to  the  U  operator  to  calculate  the  perimeter  fields  and  the  asso- 
ciated derivatives  on  the  objects  perimeter.  This  looping  may  be  repeated  until  a  steady 
state  condition,  at  point  3,  is  reached.  The  existence  of  a  steady  state  condition  assumes 
stability.  Stability  for  physical  systems  should  not  be  a  problem.  However,  when 
mathematically  modeled,  instabilities  may  result.  The  error  magnification  or  condition 
number  of  the  system  must  also  be  seriously  considered  if  this  iterative  looping  process 
is  to  be  used.    The  alternative  approach  is  to  form  an  equivalent  system,  where: 


with 


equivalent  operator  =  U »  [I  —  T *  L~\  \ 


/=  Identitv  Matrix 


and 

r»  U  =  Combined  Effects  of  the  T  and  U  operators  . 
Either  approach  is  viable  since. 

2  1 

^tocal  =  ^incideni  +  T»  U»  ij^ incident  +  {T •  11    •  ^incident  +    •    •    •     =  [/  -  T*  (.1      ^incident 

This  becomes  more  obvious  when  iA„„^„,  is  factored  from  the  equality  leaving, 

1  +  r.L^-f  (r.Lf +  .  .  . 
Placing  this  in  closed  form, 


oo 


«=0 


Finding  the  equivalent  operator  will  require  a  matrix  inversion  and  for  ver\'  large  prob- 
lems this  may  lead  to  excessive  computation  times.  The  matrix  inversion  technique  will 
be  investigated.  The  fields  may  then  be  extended  to  any  point  in  space  using  a  far  field 
Green's  function  surface  contour  integral.  This  far  field  pattern  is  available  at  point  5. 
The  incident  field  is  usually  produced  by  a  plane  wave  generator.  This  provides  the 
boundar\'  conditions  on  the  offset  contour.  This  contour  is  called  "offset"  since  it  is 
approximately  the  same  "shape"  as  the  objects  perimeter  but  is  slightly  larger.  The  dis- 
tance between  the  perimeter  and  this  contour  is  called  the  offset  distance  and  will  be 
discussed  in  Chapter  III.  The  boundary'  conditions  may  be  any  desired  field  or  wave  that 
satisfies  Maxwell's  equations.  These  waves  may  arrive  from  any  direction  and  be  of  any 
magnitude.  The  boundarv'  conditions  may  also  be  a  composite  of  any  number  of  waves 
since  superposition  does  apply  to  these  systems.  A  user  provided  subroutine  is  necessarv' 
if  conditions  other  than  a  single  plane  wave,  cylindrical  mode  (of  arbitrary  mode  num- 
ber) or  individual  input  boundary'  condition  is  desired. 

B.     MAXWELL'S  EQUATIONS 

Maxwell's  equations  can  be  written  using  our  previous  normalizations  as, 

Vx£  =  /^,77  (2.3) 

and 

V  X  77  =  L,E.  (2.4) 

[Ref  4  ]     Let  0^  =  ^— ,  with  similar  definitions  for  Dy  and  D^-    Equations  2.3  and  2.4 
can  be  further  expanded  into  the  D^,  Dy  and  D^  components  such  that. 

^,H,=  DyE,  (2.5) 

ix.Hy  =  -D^E,  (2.6) 

^l,H,  =  D^Ey  -  DyE^  (2.7) 

erE,  =  DyH,  (2.8) 

c,Ey  =  -D^H,  (2.9) 

z,E,  =  DxHy-DyH^.  (2.10) 


For  the  TM  case,  with  propagation  in  the  z  direction, 

DyE, 


//,.= 


■"■  a. 


and 


-D^E, 


Note  that  the  H^  field  =  0.   These  two  equations  can  be  combined  to  form, 

J7  =  ^V^,  X  I.  (2.11) 

Similarly  for  the  TE  case,  with  propagation  in  the  z  direction. 

£  =  J-V«/y2  X  z.  (2.12) 

Substituting  equations  2.5  and  2.6  into  equation  2.10  yields, 

Substituting  equation  2.1  into  equation  2.13  yields, 

V/',  +  /)^|^^,)  +  Z),(^^,)  =  0.  (2.14) 

Equation  2.14  can  be  further  simplified  to, 

^•[-^^'/'i]  +  ^>i=0.  (2.15) 

Similarly,  equations  2.2,  2.8,  2.9  may  be  substituted  into  equation  2.7.   This  yields, 

V.[-^V.A2]+Mr'/'2  =  0.  (2.16) 

Equations  2.15  and  2.16  are  TM  and  TE  duals.  These  two  differential  equations  describe 
the  potentials  inside  the  object  of  interest.    Defining  —  =  a  and  £,  =  /)  for  the  TM  case 


and  —  =  a  and  ^ir  =  P  for  the  TE  case  and  substituting  these  new  definitions  into 
equations  2.15  and  2.16  yields  one  differential  equation, 

V»l7.Vip^  +  Pil/  =  0.  (2.17) 

C.     VARIATIONAL  EQUIVALENCE  TO  THE  DIFFERENTIAL  EQUATION 

The  Euler-Lagrange  variational  formulation  is  based  on  the  stationarity  of  a  func- 
tional, of  the  function  \p  and  its  first  derivatives.   [Ref  4  ] 

/  =  J   I  F{X,  Y,  xjj,  Vi/^)  dXdY  (2.18) 

inside  S 

It  can  be  shown  that  the  first  variation  of  the  functional  is  zero,  61  =  0  ,  if  the 
Lagrangian.  F,  satisfies  the  Euler  -  Lagrange  equation: 

'^      ^+-TTr(   ^,i^ ,A-4^  =  0.  (2.19) 


cX  \  c(Z)^v'A)  /     -^y  \  c(Z)y(A)  J     cv> 

The  problem  thus  becomes  to  find  the  F,  which  when  substituted  into  equation  2.19, 
yields  the  original  differential  equation.   The  Lagrangian, 

can  be  simplified  to, 

/■=aVi//.V»//-/?tAl  (2.20) 

When  equation  2.20  is  substituted  into  equation  2.19, 

dF 


cF 

e[Dy4) 

cF 


-=2aZ)ytA 

=  2/?.//. 


Therefore, 


-^  (2aZ)^-.A)  +  jY  (la.DyiP)  +  2/?.A  =  0 


or, 


which  simplifies  to, 


DxluD^i/l  +  Z)y[aZ)yiA]  +  j?iA  =  0 


V  .  [aV(A]  + /?iA  =  0. 


(2.21) 


Therefore,  the  general  functional  has  been  found  since  equation  2.21  and  equation  2.17 
are  identical.    Substituting  the  a  and  /?  definitions  into  equation  2.20  yields. 


F,=^V^,.V.A,-£,0? 


(2.22) 


and 


F2  =  —  V\p2*^^h-^^r^l- 


(2.23) 


Equations  2.22  and  2.23  are  integrated  over  the  interior  to  S  with  known  boundary' 
conditions  for  either  1//,  or  (//^  on  S.  To  physically  interpret  the  variational  formulation, 
it  is  noted  that. 


and 


ViP,=^,{Hj'-HyX} 


A  A, 


WiP,  =  c,{ExY-EyX} 


Therefore, 


/i  = 


tirH  *  H  -  E,E  *  E  dX  dY 


and 


Substituting  Maxwell's  equations  into  /,  gives, 


A  = 


(V  X  E)*H-(y  X  H)*EdXdY 


s 


/,  = 


V.(£  X  H)dXdY 


s 


/,  = 


E  X  H  »ndl. 


Note  that  E  x  //  is  the  complex  oscillatory  Poynting  vector.  It  is  diflerent  from  the 
usual  Poynting  vector  oi'  E  x  H'.  Thus,  both  functionals,  /j  and  ^  are  proportional  to 
the  complex  phasors  for  oscillatory  power.  The  oscillatory  power  is  the  phasor  repres- 
enting the  excess  of  instantaneous  radiated  power  minus  the  average  radiated  power. 

D.     FINITE  ELEMENT  BOUNDARY  VALUE  SOLUTION 

With  the  discussion  of  the  variational  mathematics  completed,  a  simple  rectangular 
boundary  value  problem  will  be  discussed  in  detail.  The  rectangular  geometry  allows  for 
an  easier  formulation  but  in  no  way  limits  the  solution  from  being  extended  to  more 
comphcated  object  geometries.    Figure  3  on  page  11  shows  the  region  of  concern. 
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0 


Figure  3.      Rectangular  Object 


q      X 


Consider  spanning  the  rectangular  region  by  a  triangular  mesh,  as  shown  in 
Figure  4  on  page  12.  Both  a  and  P  may  be  functions  of  position  but  may  not  vary 
within  an  individual  triangular  element.  Given  a  fine  enough  mesh  structure,  a  smooth 
transition  in  material  properties  may  be  approximated.  For  notational  simplicity  this 
positional  dependance  will  not  be  carried  forward.  The  variational  approach  will  yield 
the  solution  to  the  boundary  value  problem  by  finding  the  ip{x,y)  which  gives  the  sta- 
tionary value  of, 


na 


1  = 


rb 


(aViA  .  ViA  -  /?iA^)  dx  dy 


*'0    •'0 


which  is  constrained  on  the  boundary  by  the  previously  specified  boundary  conditions. 
[Ref  5  ] 
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N+1 
N 


J   =   0 

i    -   0  1 

Figure  4.      Rectangular  Mesh 


M         M+1 


1  he  values  of  fp  at  the  interior  nodes  become  the  discrctized  unknowns: 


where 


^..  =  ^(a;  }})     i=\...M  and  J  =\  ...N 


A)  =  i»AX 


and 


Yj=j.AY 


for  the  uniform  mesh  structure  with  AX  = 


{M+D 

A/4-1  A'+l 


and  AY  = 


(N+l) 


.  Approximating, 


^{x,y)==  Yj  Yj'^iJ^iM^y'^ 


/=0    7=0 


which  includes  the  known  boundary  nodal  values  (/^o.y  and  ^m^^j  for  j  =  0  and  N+l, 
and  i/',  0  and  ^i^s^x  for  i  =  0  and  M  +  1,  linear  pyramidal  basis  functions,  u,j{x,y),  which 
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have  unit  value  at  the  (i,j)  node  and  zero  value  at  all  surrounding  nodes.    Figure  5(a)  is 
a  top  view  of  the  pyramidal  basis  function,  while  Figure  5(b)  is  a  perspective  view. 

1.0 


a.  b. 

Figure  5.      Pyramidal  Basis  Function,  (a)  top  view,  (b)  perspective  view 

Ihc  functional  1  will  thus  be  a  discrete  function  of  each  of  the  nodal  values  of  ^, 

/=  ^("Aco.  "Ao,  1 'A/,;. -.'/'w+i.^+i)- 

The  approximate  discrete  solution  will  be  found  by  the  system, 


81 


^^m.n 


=  0,     for    m  =  1  .  . .  M  and  «  =   1  .  . .  A'. 


Now, 


ei 


S^m.n 


=  2 


ra    (*b 


■'0     ''O 


a  ..         »Vi//-/?  ^/      lA  ]dxdy  =  0 


^^m,n 


^^m.n 
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where, 


a;+i  a'+i 


V^{x,y)=Y,  Yj'^ij'^^^iM^y)- 


i=0    >0 


The  gradient  of  the  basis  function  is, 


^"^,X^.j)  = 


and  the  basis  function  is, 


^m,ni-'^^y)  =  -j:;, 


exit 


(^^m,n 


Therefore,  the  system  of  equations  to  solve  becomes, 


(aVz/„,^  „  .  Vi//  -  fh(„^,^\lj)  dx  dy  =  0,  for  m  =  1  ...  A/  and  n=  \  .  .  .  N. 


Substituting      for      Vi//  and  \p      in      terms      of     the      nodal      values      of     ij/      for 
m  =  I  .  .  .  M  and  n  =  \  .  .  .  X  gives. 


M+\  .V+1 


{c^.^u^^  „  .  Vuij  -  (iu^^  „uij)  dx  dy  =  0 


i=i  y=o 


where, 


i^-^^m,  n  •  V«/,y  -  /?w^,  „W;-.y)  t/x  dy  =  F{{m,  «),(/, y)} . 

Regrouping  this  to  put  the  known  nodal  values  on  the  right  hand  side  gives  for 
m  =  \  .  .  .  M  and  /;  =  1  .  . .  A"  (interior  nodes). 


A/      .V 

IE- 

/=1  7=1 


il^yFlim,  «),(/,;)]  =  -      ^      Y.'^yFlim,  n),(/,y)]. 


Boundary 

Sodcs  Only 

for  (ij) 
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By   renumbering   the   nodes    using   a    single   index   for   the   unknown   nodal   values, 
k  =  N{i  -  \)+J  and  /  =  Xini  -  \}  +  n   and, 


M'\ 


yF{i,k)ii^j,=-yni,k')^,, 


k=\ 


k' 


The  functional  F{1,  k)  =  0  if  nodes  /  and  k  are  not  both  associated  with  at  least  one 
common  triangular  element.  Therefore,  Vu^  •  Vw^  and  UiU^  will  be  zero  except  in  triangles 
where  /  and  k  both  appear  as  nodes.  For  the  example  mesh  structure  of  Figure  3  on 
page  11,  this  produces  a  banded  matrix.  F.  This  sparse  matrix  can  be  easily  displayed 
by  first  defining  column  vectors  of  the  unknown  nodal  values  in  the  mesh 

"^1=1^1,1^  '/'/,2'--M  "A/.  ,v]^- 

The  F  -  matrix  elements  F  [(m,  «),(/,y)]ire  zero  unless  the  (m,n)  node  shares  at  least  one 
element  with  the  (i.j)  node.  Thus,  the  node  values  in  \p,  will  be  coupled  only  to 
(//,_,,  t//,.,  and  any  associated  boundary  nodes.   This  is  written  as, 

where  A„  B,  and  C,  are  A'  x  A'  complex  arrays  and  P,  is  a  A'  x  \\  array  where  A'^  is  the 
number  of  boundary  nodes  associated  with  the  i-th  column.   This  appears  as, 


B^ 

c, 

^2 

Bi 

Q 

^3 

B3 

C3 

A 


:V/-1 


- 

- 

- 

^^ 

^1 

'''', 

'^'2 

P2 

'*'*: 

0 

M'3 

P3 

'•'«, 

• 

• 

• 

• 

• 

• 

• 

• 

• 

V, 

^M-\ 

^A/-, 

Pm-1 

'''b.u-i 

Am 

Bm 

^M 

Pm 

">'«« 

. 

_ 

- 

. 

If  we  denote  t/^o  as  the  initial  boundary  values  and  i/'w-i  ^^  the  final  boundary  values 
then  i//£^  becomes  only  the  boundary  conditions  on  the  top  and  bottom  at  j  i^  0  and 
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N  +  1.  Note  that  the  above  system  is  tri-block  in  nature  and  has  a  large  number  of  zero 
valued  elements.  The  zero  valued  elements  were  omitted  for  clarity.  Each  element  is 
actually  a  matrix  and  therefore  it  is  evident  how  sparse  this  system  is.  Each  of  the 
A„  B,.  C,  and  Pj  matrices  are  equally  sparse,  however,  a  global  symmeti^  does  not  ap- 
pear. 

E.     EVALUATION  OF  THE  F  -  MATRIX  CONTRIBUTIONS 

Given  an  arbitrary  element  as  shown  in  Figure  6,  the  potential,  i//,  can  be  linearly 
approximated  by  [Ref  5  ], 


«/'(-v,>)  =  (.r,j-,  i).[r] 


3 
k=\ 


^3 


Figure  6.      Mesh  Element 

where  t//^  =  (//(-V;^,^;)  is  the  nodal  value  at  the  k-th  node,  [T]  =  3  x  3  Transform  Array 
and 
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Ui^{x,y)  =  Linear  basis  functions  for  the  k-th  node 


=  ix,y,l) 


.k 


T 


2,k 


—  '^\,k^  +  ^2,  ^>'  +  ^3,  k- 


Note  that, 


"/{(•^w  J'w)  — 


1    ,  for  k  =  m 

0    ,  for  k  i=^  m 


It  can  be  shown  that. 


[r]  = 


2A 


Lvi-y^)  (j'3-Ji)  Cvi-j-2) 

(-Y3--^2)  (-^l--^3)  (-^2-^]) 

{x.Ji  -  X^V2)         {X^\\  -  X^Vi)         (.V,J-2  -  X,J\] 


where  M  |  =  triangle  area,  and 


2A  =  det 


'^1 

y\ 

r 

-^2 

y'l 

1 

-^3 

>"3 

1 

2 A  =  ix^V2  +  x^V]  +  x^yj)  "  (•^3>'2  +  ^i>3  +  -^2^1) 


Furthermore,  the  linear  basis  functions,  Ui,{x,y),  can  be  interpreted  as  relative  areas  of 
the  triangle  shown  in  Figure  7  on  page  18. 
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Figure  7.      Basis  Functions  Interpretation 


A  =  A^  +  A2  +  A2 


is  constant  as  (jc,>0  varies  and, 


«^a(-^,>')  = 


Akix,}') 


Within  a  given  triangular  element,  the  evaluation  for  k  =  1,  2,  3  and  /  =  1,  2,  3  in  the 
element  is  of  interest  and, 


inside 
triangle 


{aS^Ui .  Vuj^  -  pMUi,)  dx  dy. 


The  assumption  is  made  that  a  and  /?  will  be  approximated  as  constants  within  the  tri- 
angle. These  material  constants  can,  however,  vary  from  element  to  element.  Taking 
the  gradient  terms  first, 
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Therefore, 


1 


o.yui  •  Vuifdxdy 


triangle 

where  A  is  the  area  of  the  q"'  element.   Next  it  can  be  derived  that, 


"        [  -^\A\     ,  fork  =  l] 

triangle 

More  generally,  with  each  u,  raised  to  an  integer  power.  n„ 


u"^U2'U2^dxdy  =  2\A\ 


^h'"2-^b- 


tnansle 


(«i  +  th  +  n^  +  2)! 


The  final  result  is. 


Fail,  k)  = 


{y.\'ui  •  Vui^  —  jJuiUi^)dxdy 


triangle 

q 


=  \A\^,.{Tl,+  Tl,)-^py     k  =  L 

F.     GREEN'S  FUNCTION  CONTOUR  INTEGRAL 

The  scattered  fields,  ^^  from  an  arbitrary  object  in  a  vacuum  satisfying  Helmholtz's 
equation  (see  Figure  8  on  page  20  ), 


vV  +  ^^»A  =  o. 


are  [Ref  6  ], 
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•PA') 


G(HP)  "'' 


on 


on 


contour 


Figure  8.      Green's  Function  Integration 


where  the  Green's  function  is, 


G(r|F')  =  4<(*ol^--r'l) 


and 


dn 


=  n  *V\Ij     on  the  contour 


and 


dG       A    A  y^o   ,,(!)( ,   I-      -,  I  \ 
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The  Hankel  functions  of  the  second  kind  of  order  zero  and  one,  H^^^  and  //|^'  ,  will 
present  a  problem  for  the  numeric  integration  discussed  in  Chapter  V.  The  imaginary 
portion  of  these  functions  rapidly  approaches  negative  infmity  as  the  argument  ap- 
proaches  zero.     The  -r: —  is  obtained  bv  a  finite  difference  method  usins  the  field 

CN 

boundary  conditions  on  surface  B  (boundary  conditions)  and  the  calculated  field  condi- 
tions on  surface  P  (object  perimeter).   This  results  in, 

Oy/  'r  boundary        Y  perimeter 

cn   ^       offset  distance 

It  will  be  shown  that  to  maximize  the  accuracy  of  the  Green's  function  the  offset 

c\l/ 
distance  should  be  made  as  laree  as  possible.    This,  however,  causes  the  -r; —  to  be  m- 

cn 

accurate.  Thus,  an  optimal  condition  must  be  found  that  maximizes  the  accuracy  of  the 
entire  numeric  integration.  Such  a  condition  does  not  maximize  the  accuracy  of  any  one 
of  the  contributing  parts  to  the  Green's  function  integrand. 

G.     FAR-FIELD  EVALUATION 

When  the  Green's  function  integral  discussed  above  is  used  for  far  field  calculations 
several  simplifying  relationships  develop.  These  simplifications  require  a  less  demanding 
numerical  integration.   To  be  in  the  far  field  region  three  conditions  must  exist, 

|F|  >Z).         |/"l  >/n       and        |F|  >  ^^ 


^■0 


where  /q  is  the  free  space  wavelength  and  D  is  the  maximum  dimension  of  the  object. 
[Ref  7  ]  As  X  approaches  infmity. 


");■ 


//«>(.v)^V^.-^' 


and 


"); 


Hf\.)^^^je-'\ 


This  requires, 


^('^|'"")47i|^^"*'' 


and 
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cG 


en 


=  —  n  •  K  —r 


In  the  far  field 


1  /  1 


Therefore, 


and 


V  R 


4    V   kKqR 
yj  -jr  and  h  •  R^  h  *  r.   Thus, 


—  .-^•^°^. 


/^f-\-'\  ^         /      '^''  — y^n'' -/'^n'^  cos  ( 


('^i'^')  =  -4y^^"*("-'V' 


With  these  new  definitions   substituted   into   the   original  Green's  function  integral 
equation, 


^scaiiercJ*') 


Cl//  A  A         ,_ 


^Jk,r' cos  e^^, 


Note  that  this  equation  is  partitioned  into  a  distance  dependent  term  and  a  theta  depend 
term.   The  theta  dependent  term  may  be  defined  as 


on 


^■V'cos^^^,_ 


The  two  dimensional  bistatic  radar  cross  section  (RCS)  per  unit  length  of  the  cylindrical 
structure  may  now  be  defined  as, 


2nrP'       ,.       2nr\\l/'\ 


5|2 


RCS  =  a((b\  S')  =  lim  — — : —  =  lim 

V^    '  ^    /       r-^co       ni'^c  r-oc 


oo      pine  r_oo         I  ^L'  I  2 
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a  =  lim  2nf 


Snhr         4kQ 


2n 


where   the  wavenumber,   A'o  =  -^^—  and  the  incident   field  is  assumed  to   be   of  unit 
magnitude,  1 1//'|  =  1.0. 
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III.     MESH  GENERATION 

A.  INTRODUCTION 

A  display  or  plot  of  the  computer  generated  mesh  structure  is  not  required  for  the 
problem  solution,  however,  it  provides  an  immediate  visual  confirmation  that  the  in- 
tended problem  geometr\'  has  been  entered  correctly  into  the  computer.  A  large  amount 
of  initialization  data  is  required  for  even  the  most  rudimentarv'  problem.  For  this  reason, 
the  input  data  is  provided  to  the  mesh  generation  program  via  a  data  file.  Modifications 
are  possible  at  a  later  time  with  minimum  efibrt. 

B.  INPUT  DATA 

The  input  data  file  is  called  IXPUT.DAT  and  contains  25  input  fields.  Of  these 
fields,  14  relate  directly  to  the  generation  or  display  of  the  mesh  structure.  Only  the 
object  surface  coordinates  require  adherence  to  a  specific  format.  All  other  data  need 
only  be  of  the  correct  type  (i.e..  character,  integer  or  real).  A  brief  description  of  each 
field  is  provided  below. 

•  Field  1  is  a  label  of  no  more  than  12  characters.  This  label  is  for  the  plot  and  input 
data  file. 

•  Field  2  is  a  character  flag  that  specifies  the  input  coordinate  system.  If  set  to  "R" 
or  "r".  the  rectangular  system  is  used.  If  set  to  "P"  or  "p".  the  polar  system  is  used. 
If  a  "P"."p"  ."R"  or  "r"  is  not  detected,  an  error  is  returned  to  the  display. 

•  Field  3  is  a  character  flag  that  if  set  to  "V  or  "i"  will  cause  several  intermediate 
values  to  be  stored  to  disk  during  the  Finite  Element  Boundary  Value  (FEBV) 
program  execution.   This  option  was  used  during  the  debug  process. 

•  Field  4  is  a  character  flag  that  if  set  to  "D"  or  "d"  will  create  a  DISSPLA 
FORTR.AN  program  capable  of  replicating  the  input  object,  in  wavenumber  nor- 
malized coordinates,  on  mainframe  computers  having  a  DISSPLA  graphics  pack- 
age. DISSPLA  is  a  subroutine-based  language.  The  generated  program,  called 
DISSPLA. FOR,  is  a  compilation  of  four  subroutines  calls  per  element.  This  file 
can  get  very  large  for  dense  mesh  structures. 

•  Field  5  is  a  character  flag  that  if  set  to  "L"  or  "u"  will  cause  a  uniform  material  to 
be  assumed.  In  the  uniform  case,  no  material  interface  exists.  This  option  was 
only  used  to  verify  the  Finite  Element  solution  accuracy. 

•  Field  6  is  a  character  flag  that  if  set  to  "M"  or  "m"  will  cause  only  the  mesh  to  be 
generated.  This  is  very  useful  when  first  starting  a  problem  and  the  optimal  mesh 
structure  has  not  been  determined. 

•  Field  7  is  a  real  number  specifying  the  desired  mesh  resolution  in  wavelengths.  The 
mesh  resolution  determines  the  dimension  of  the  mesh  elements. 
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Field  8  is  a  real  number  specifying  the  distance,  in  wavelengths,  between  the  object 
perimeter  and  the  ofTset  boundar}'  contour. 

Field  9  is  a  real  number  that  specifies  a  multipHcative  scaling  factor  for  the  nu- 
merical integration  stepping  function  discussed  in  Chapter  V. 

Field  10  is  a  bias  term  that  can  be  used  to  shift  the  numerical  integration  stepping 
function.   This  term  is  used  for  distances  less  than  1.0. 

Field  1 1  is  a  bias  term  that  can  be  used  to  shift  the  numerical  integration  stepping 
function.   This  term  is  used  for  distances  greater  than  1.0. 

Field  12  is  a  real  number  specifying  the  maximum  distance  beyond  which  no  further 
contribution  to  the  Green's  Function  Integral  is  made.  If  this  feature  is  not  de- 
sired, this  term  should  be  made  larger  than  the  objects  maximum  dimension  plus 
twice  the  ofiset  distance. 

Field  13  is  an  integer  specifying  the  number  of  input  data  points. 

Field  14  is  an  integer  specifying  the  angular  resolution,  in  degrees,  desired  for  the 
final  radar  cross  section  calculation. 

Field  15  is  an  integer  specifying  the  mesh  generation  technique. 

Field  16  is  an  integer  specifying  the  perimeter  node  from  which  the  bisection  seg- 
ment originates.   This  node  is  called  the  "start  node". 

Field  17  is  an  integer  specifying  the  perimeter  node  on  which  the  bisection  segment 
terminates.    This  node  is  called  the  "stop  node". 

Field  IS  is  a  pair  of  real  numbers  (on  two  hnes)  specifying  the  x  and  y  coordinates 
by  which  the  object  will  be  displaced. 

Field  19  is  a  real  number  specifying,  in  wavelengths,  the  desired  distance  between 
the  origin  and  the  first  input  data  point.  Fields  18  and  19  when  used  together,  al- 
low an  object  to  be  placed  at  any  position  and  scaled  to  any  size. 

Field  20  is  a  pair  of  real  numbers  (on  two  lines)  specifying  the  real  and  imaginar>' 
parts  of-—  for  the  TM  case  or  7-  for  the  TE  case. 

Field  21  is  a  pair  of  real  numbers  (on  two  hnes)  specifying  the  real  and  imaginary 
parts  of  £,  for  the  TM  case  or  n,  for  the  TE  case. 

Field  22  is  a  character  fiag  that  if  set  to  "P"  or  "p"  enables  a  plane  wave  generator. 
The  plane  wave  is  propagating  down  the  y  axis,  and  generates  an  E^e'''^  condition 
on  the  offset  boundary  contour.  If  the  flag  is  set  to  "C"  or  "c",  a  cylindrical  mode 
generator  is  enabled.  This  generates  an  E^  cos  n<i)  condition  on  the  offset  boundary 
contour.  If  a  "P"  ,  "p",  "C"  or  "c"  is  not  detected,  then  manually  input  boundary- 
conditions  must  follow,  and  fields  23  and  24  are  not  used. 

Field  23  is  a  real  number  specifying  the  wave  amplitude. 

Field  24a  is  a  real  number  specifying  the  wave  frequency  (in  Hz). 

Field  24b  (only  for  cylindrical  case)  is  an  integer  specifying  the  mode  number,  n. 

Field  25  is  the  object  perimeter  data,  in  either  polar  or  rectangular  form. 


25 


An  example  of  an  input  data  file  for  a  homogeneous  circular  cylinder  is  provided  in 
Appendix  A.  The  majority  of  the  input  data  is  echoed  to  the  computer  display  and  a 
system  "pause"  is  initiated  to  allow  for  user  inspection.  The  program  may  be  aborted 
or  continued  at  this  time.  The  initial  object  dimensionalization  has  already  occurred, 
and  the  number  of  unknowns  and  the  maximum  unknown  width  is  displayed.  These 
factors  give  an  excellent  indication  of  the  expected  run  time  for  the  FEBV  routines.  For 
example,  a  problem  with  512  unknowns  and  a  maximum  unknown  width  of  31  took  806 
seconds  to  execute  while  a  run  with  8  unknowns  and  a  maximum  unknown  width  of  3 
took  only  10  seconds  to  execute.  These  times  are  for  a  Intel  80386  based  personal 
computer  with  an  Intel  80287  co-processor  chip  calculating  the  fields  for  a  circular  cyl- 
inder.  The  FEBV  routines  are  the  next  code  block  to  execute  after  the  pause  is  cleared. 

C.     MESH  GENERATION  PROGRAM 

The  mesh  generation  program  consists  of  seven  subroutines.  These  subroutines  are 
an  integral  part  of  the  finite  element  program  and,  therefore,  were  not  separated.  These 
routines  are  discussed  below. 
1.     lO  (Input/Output) 

This  subroutine  reads  the  information  contained  in  the  INPUT.DAT  file  dis- 
cussed earlier.  A  two  dimensional  object  can  be  described  in  any  number  of  ways, 
however,  for  simphcity,  the  polar  and  rectangular  coordinate  systems  are  used.  In  either 
case,  the  initial  assumption  is  that  all  data  points  are  referenced  to  a  local  origin.  This 
local  origin  can  be  offset  by  any  desired  amount  using  field  18.  This  offset  is  independ- 
ent of  the  entered  data  points  or  any  size  scaling  provided  by  field  19.  A  plot  label  file, 
named  TEXT.LBL,  is  created  and  the  initial  object  perimeter  (coded  for  a  display  in 
blue)  is  written  to  the  output  file,  PLT.DAT.  These  data  points  describe  the  perimeter 
of  the  object.  This  will  later  prove  helpful  in  determining  the  conformity  of  the  gener- 
ated mesh  to  the  input  perimeter.  All  subsequent  screen  writes  are  coded  for  a  display 
in  green.  Two  example  objects  are  shown  in  Figure  9  on  page  27.  These  objects  will 
be  used  throughout  this  chapter.  The  circular  cylinder  was  generated  by  a  separate 
computer  program.  The  "horseshoe"  shaped  object  was  manually  input.  Graph  paper 
was  used  to  determine  the  x  and  y  coordinates  of  the  28  unequally  spaced  perimeter 
points. 
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OBJECT    PERIMETERS 

Figure  9.      Typical  Objects 

2.  Rotate 

This  subroutine  reorders  the  input  data  points  to  allow  for  any  desired  bisection 
segment  start  and,  or  stop  node.  This  subroutine  is  only  used  if  manual  selection  of  the 
bisection  segment  start  and; or  stop  nodes  is  requested. 

3.  Bound  (Boundary) 

This  subroutine  sub-divides  the  object  perimeter  based  on  the  mesh  resolution 
specified  in  Field  7.  The  mesh  resolution  is  the  approximate  length,  specified  in  wave- 
lengths, that  the  user  desires  the  perimeter  to  be  divided  into.  The  division  of  the  per- 
imeter is  based  on  linear  interpolation  between  input  data  points.  A  new  perimeter  node 
is  placed  at  each  of  the  sub-division  points.  Additional  input  data  points  are  necessary 
in  areas  of  rapid  change  to  allow  for  a  correct  object  perimeter  representation.  Ihe  ob- 
ject bisection  extends  from  the  "start  node"  to  the  "stop  node".  These  nodes  are  user 
specified  in  Fields  16  and  17  and,  in  general,  divide  the  object  in  half  Ihe  bisection 
should  be  arranged  such  that  the  object  width,  perpendicular  to  the  bisection  segment, 
is  minimized.  Thus,  a  long  slender  object  should  be  oriented  for  a  major  axis  bisection. 
Whether  the  major  axis  is  oriented  vertically  or  horizontally  is  of  no  importance.  This 
is  demonstrated  in  Figure  10  on  page  28  (right  side). 
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BISECTION    SEGMENT 


Figure  10. 


•—OBJECT    PERIMETER-' 

Bisection  Segment  and  Row  Arrangement 


For  more  complex  objects,  the  bisection  is  not  a  straight  line,  but  rather  a  series 
of  line  segments  that  approximate  a  curve.  The  number  of  perimeter  nodes  on  the  left 
side  of  the  bisection  segment  must  equal  the  number  of  perimeter  nodes  on  the  right  side 
of  the  bisection  segment.  The  bisection  segment  also  contains  this  same  number  of 
nodes.  Thus,  the  program  must  adjust  the  requested  mesh  resolution  to  ensure  proper 
nodal  spacing.  This  allows  for  a  piecewise  continuous  segment  to  cross  the  object. 
4.     Normal 

This  subroutine  calculates  the  unit  normals  to  the  object  perimeter  at  each 
newly  established  perimeter  node.  The  normal  is  a  perpendicular  constructed  to  the 
chord  connecting  the  two  nodes  adjacent  to  the  node  for  which  the  calculation  is  oc- 
curring. This  perpendicular  originates  at  the  current  node  and  is  of  unit  length.  See 
Figure  1 1  on  page  29. 
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UNIT    NORMAL 
CURRENT    NODE 


•—ADJACENT    NODES 


Figure  11.      Unit  Normal  Calculation 

5.     Nodset  (Node  Set) 

This  subroutine  uses  a  two  sweep  technique  to  compute  the  number  of  nodes 
on  each  nodal  row.  The  rows  are  labeled  I  =  1,  2,  3,  ...,  I  maximum  (I MX)  consec- 
utively from  the  top  of  the  object  to  the  bottom.  Segments  normal  to  the  perimeter 
surface,  calculated  in  NORMAL,  connect  the  perimeter  nodes  to  the  offset  boundary 
contour.  Two  such  completed  normal  segments,  one  on  each  end,  complete  a  nodal  row. 
When  all  rows  are  complete,  a  second  contour  has  been  established  that  approximates 
the  object's  perimeter.  This  contour  is  displaced  from  the  perimeter  by  the  distance 
specified  in  Field  8.  This  contour  provides  the  boundary  conditions  for  the  FEBV  rou- 
tines.  The  optimal  selection  of  this  offset  distance  will  be  a  topic  of  Chapter  IV. 

With  the  object  divided  into  rows,  each  row  can  be  further  divided  into  equally 
spaced  nodes.  Based  on  the  requested  resolution  and  whether  the  objects  dimensions 
are  expanding  or  contracting,  the  nodal  spacing  is  adjusted  to  keep  the  elements 
approximately  the  same  size.  Two  adjacent  nodal  rows  form  an  elemental  row.  These 
rows  are  given  the  same  label,  I  =  1,  2,  3,  ...,  I  MX  as  the  upper  nodal  row.  The  nodes 
on  an  elemental  row  are  numbered  consecutively,  starting  with  the  left-most  node.  1  his 
node  is  always  an  offset  boundary  contour  node.     When  the  upper  nodal  row  is 
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numbered,  the  process  is  continued  for  the  lower  nodal  row.  An  example  of  this  element 
row  numbering  scheme  is  shown  in  Figure  1 2  on  page  30. 


Figure  12.      Element  Ron  Numbering  Scheme 

A  mesh  orientation  attribute  is  set  for  the  left  and  right  portion  of  each  element  row. 
This  attribute  determines  if  the  object  has  a  clockwise  or  counter-clockwise  orientation. 
See  Figure  13  on  page  31.  Switching  between  the  four  available  mesh  orientations  al- 
lows for  a  mesh  that  more  accurately  conforms  to  the  input  object  perimeter  without 
resulting  in  a  disproportionate  mesh  structure. 
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LEFT    SIDE 


RIGHT    SIDE 


MOR    =    1     MOR    =    0 


MOR    =    D     MOR    =    1 


Figure  13.      Mesh  Orientation  Attributes 

There  is  an  additional  requirement  that  the  number  of  nodes  on  a  given  left  or  right  half 
row  must  be  only  one  more,  equal  to,  or  one  less  than  the  number  of  nodes  on  the  ad- 
jacent left  or  right  half  row.  Thus,  a  two  sweep  process  is  utilized.  This  process  ensures 
that  this  requirement  is  met  and  that  the  first  and  last  row  have  only  two  nodes.  The 
second  and  second  from  last  row  (if  present)  must  have  three  nodes.  It  is  possible,  as 
shown  in  Figure  14  on  page  32,  to  generate  a  mesh  that  has  only  three  rows.  This  object 
uas  input  as  a  circular  cylinder.  It  is  evident  that,  due  to  a  small  number  of  elements, 
the  generated  mesh  more  accurately  resembles  a  square.  This  will  be  a  problem  for  cal- 
culating the  scattered  fields,  however,  the  internal  fields  can  be  accurately  approximated. 
In  this  special  case,  the  second  and  second  from  last  rows  are  the  same.  It  will  be  shown, 
in  Chapter  V,  that  since  this  mesh  does  not  closely  approximate  the  input  objects  per- 
imeter the  resulting  Green's  function  contour  integrals  and  subsequent  far  field  calcu- 
lations have  reduced  accuracy. 

The  nodes  of  the  nodal  rows  form  the  vertices  of  triangular  elements.  An  ele- 
ment, as  defined  in  Chapter  II,  has  three  vertices,  each  assigned  a  unique  number  (I,  2 
or  3).  The  ordering  of  these  vertices  depends  on  the  mesh  orientation  attribute.  Each 
element  also  has  an  associated  relative  dielectric  constant  (£,)  and  relative  permeability 
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ROW    1 


ROW    2 


ROW    3 


Figure   14.      Three  Ron  Cylinder 


6.     Sorter 

This  subroutine  generates  a  complete  mesh  row  and  all  the  element/node  inter- 
connection relationships.  With  the  nodal  structure  in  place,  individual  nodes  can  be 
assigned  to  individual  elements  within  the  global  mesh  structure.  Each  element  in  a 
given  row  is  assigned  a  unique  local  element  number  starting  with  the  left  most  element 
(relative  to  the  bisection  segment).   See  Figure  15. 


Figure  15.      Element  Numbering  Scheme 
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It  is  vital  that  a  method  of  determining  which  nodes  form  the  vertices  of  a  given  element 
and  which  elements  have  a  vertex  attached  to  a  given  local  node.  The  nodes  that  are 
not  part  of  the  offset  contour  have  unknown  field  values.  A  single  node  can  be  con- 
nected to  as  many  as  six  elements,  only  four  of  which  can  be  in  the  current  row.  Ihis 
relationship  is  shown  in  Figure  16(a).  This  hexagonal  arrangement  of  six  elements  is 
very  common  within  the  mesh.  Along  the  bisection  segment  or  where  the  mesh  orien- 
tation attribute  changes  from  one  row  to  another,  this  pattern  is  disrupted.  An  example 
of  an  extreme  case  (the  center  of  a  circular  cylinder)  is  shown  in  Figure  16(b).  After  all 
elements  and  nodes  in  an  elemental  row  are  assigned,  an  ordered  sweep  is  conducted  of 
that  elemental  row  to  determine  which  nodes  are  connected  to  which  elements,  which 
elements  are  connected  to  which  nodes,  and  how  many  elements  a  single  node  is  con- 
nected too.   A  complete  row  has  now  been  generated.    See  Figure  17  on  page  34. 


4  ELEMENTS 
ATTACHED  TO 
THIS  NODE 


3  ELEMENTS 
ATTACHED  TO 
THIS  NODE 


Figure  16. 


T>>o  Possible  Element  Intersections,  (a)  extreme,  (b)  normal 


The  information  necessary  to  generate  the  rows  and  elements  will  be  used  again,  in 
Chapter  IV,  to  solve  the  (FEBV)  problem. 
7.     Finder 

This  subroutine  determines  the  x,  y  coordinates  of  each  node.  This  data  is  also 
necessary  to  solve  the  FFBV  program  of  Chapter  IV  and  provides  the  plotting 
coordinates  for  the  PLT.DAT  file.  1  his  process  is  repeated  for  all  rows.  Thus,  the  entire 
object  is  generated  as  the  compilation  of  elemental  rows  made  of  triangles.  See 
Figure  18  on  page  34. 
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BISECTION  SEGMENT 


OBJECT  PERIMETER 


OFFSET  CONTOUR 


Figure  17.      Typical  Row  Structure 


Figure   18.      Global  Mesh  Structure 

A  file  is  now  available  for  display  using  a  commercially  available  program  called 
"CURVE-DIGITIZER".  Any  program  that  can  accept  x,  y  coordinate  data  will 
accomplish  the  same  display  process.  Minor  changes  to  the  MESH  program  provided 
in  Appendix  B  may  be  necessary  since  several  "CURVE-DIGITIZER"  plotting  codes 
are  embedded  in  the  file  generation  code. 
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D.     OPTIMIZATION  OF  THE  MESH 

For  most  objects,  it  is  recommended  that  the  MESH  program  be  executed  in  the 
mesh  generation  only  mode  (Field  6  set  to  "M"  or  "m")  prior  to  the  execution  of  the  total 
finite  element  program.  This  will  ensure  that  the  desired  mesh  structure  is  obtained  prior 
to  solving  for  the  unknown  field  values.  The  MESH  (generation  only)  program  requires 
only  a  few  seconds  for  even  the  most  dense  mesh  structures.  It  is  readily  seen  that  as 
the  mesh  resolution  increases,  the  number  of  rows  increases  linearly  while  the  number 
of  unknowns  increases  geometrically.  Therefore,  the  mesh  calculation  times  may  not 
change  appreciably  when  the  mesh  is  made  more  dense,  however,  the  calculation  time 
for  the  finite  element  program  will  rise  geometrically.  For  this  and  other  reasons,  the 
mesh  density  should  be  kept  low.  This  will  lead  to  a  smaller  number  of  unknowns  and 
result  in  faster  program  execution  times.  Figure  19  on  page  36  plots  this  relationship 
for  a  circular  cylinder.  The  solution  for  these  unknown  field  values  will  be  a  topic  of 
Chapter  IV. 

Six  different  mesh  generation  methods  are  available  to  create  mesh  structures.  Some 
of  these  methods  were  evolutionary'  in  nature  and  provide  limited  practical  benefit. 
Methods  1  and  6  are  by  far  the  most  useful.    Each  method  is  discussed  below. 

Method  1  constructs  the  bisection  segment  by  connecting  the  first  input  data  point 
to  the  midpoint  of  the  perimeter.  This  segment  is  divided  into  equal  length  segments 
each  separated  by  a  node.  This  method  is  useful  for  ver\'  simple  objects  such  as  the 
circular  cylinder  shown  in  Figure  20  on  page  37.  This  circular  cylinder  will  be  used  in 
the  explanations  of  all  mesh  generation  methods.  These  illustrations  are  not  intended 
to  optimize  the  mesh  structure  but  rather  allow  for  easy  comparison  of  the  6  basic 
methods. 

Method  2  constructs  the  bisection  segment  as  described  in  method  1.  The  bisection 
segment  nodes  are.  however,  ordered  differently.  This  segment  is  divided  by  connecting 
line  segments  from  corresponding  nodes  on  the  left  and  right  side  perimeter.  Where 
these  line  segments  intersect  the  bisection  segment,  a  node  is  placed.  This  leads  to  un- 
equally spaced  bisection  segments.   This  can  be  seen  in  Figure  21  on  page  37. 
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Figure  19. 


BIsaellen  Nods  Numbsr 

Unknowns  and  Number  of  Rons  Versus  Mesh  Resolution 


Method  3  modifies  method  1  by  allowing  the  user  to  specify,  using  Field  17,  the  stop 
node  for  the  bisection  segment.  Ihis  method  can  be  useful  to  rapidly  adjust  around 
slightly  irregular  objects.    This  can  be  seen  in  Figure  22  on  page  38. 

Method  4  combines  methods  2  and  3  by  using  the  connected  line  segment  technique 
of  method  2  to  determine  the  node  positions  on  the  bisection  segment,  the  stop  node  of 
which  is  specified  as  in  method  3,   This  can  be  seen  in  Figure  23  on  page  39. 

Method  5  improves  upon  method  4  by  repositioning  the  unequally  spaced  bisection 
nodes.  Linearly  interpolated  positions  for  the  nodes  leads  to  equal  spacing.  This 
method  reduces  the  node  "bunching"  that  frequently  occurs  with  methods  3  and  4.  This 
can  be  seen  in  Figure  24  on  page  40. 
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Figure  20.      Method  1  Mesh  Structure  Example 


Figure  21.      Method  2  Mesh  Structure  Example 
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Figure  22.      Method  3  Mesh  Structure  Example 

Method  6  is  the  final  improvement  in  which  metiiod  5  was  modified  to  allow  for  a 
user-specified  start  node.  This  provides  the  user  with  the  ability  to  start  and  stop  the 
bisection  segment  at  any  input  node  without  having  to  rearrange  the  input  data.  Since 
method  6  contains  all  of  the  capabilities  of  the  other  five  methods,  it  is  almost  exclu- 
sively used  for  mesh  generation.  This  can  be  seen  in  figure  25  on  page  40.  There  are 
situations  that  could  be  best  served  by  one  of  the  other  methods.  For  example,  a  cir- 
cular cylinder  and  other  simple  symmetric  objects  can  be  represented  using  method  1. 
Method  2  would  work  well  with  a  square  or  diamond  shape.  Ihese  objects  would  be 
bisected  by  a  diagonal.  Method  3  would  be  most  suited  for  a  symmetric  object  with  a 
planar  material  interface.  The  more  dense  mesh  would  be  used  for  the  higher 
permittivity  material.  Method  4  would  work  well  with  a  rhombus  or  other  slightly 
asymmetric  object.  Method  5  is  almost  as  versatile  as  method  6,  and  would  work  well 
in  any  situation  where  the  start  node  is  fixed.  A  great  deal  of  planning  is  not  needed  in 
designing  most  mesh  structures  since  they  are  calculated  and  displayed  in  a  matter  of 
seconds. 
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Figure  23.      Method  4  Mesh  Structure  Example 

Iterative  selection  of  diirerent  mesh  generation  parameters  has  proven  to  be  the  best 
technique.  A  summary  of  all  mesh  generation  capabilities  is  provided  in  Table  1. 

Appendix  C  contains  a  program  called  READ. FOR.  This  program  takes  the  output 
data  file  from  the  "Curve-Digitizer"  CAD  program,  called  FINALDWG.DAT,  and  after 
receiving  the  answers  to  several  prompted  questions,  creates  a  new  INPUT. DAT  file. 
Typically,  the  CAD  program  is  used  to  generate  the  perhneter  of  the  object.  1  his  may 
be  as  a  series  of  points,  line  segments  or  a  combination  of  the  two.  The  answers  to  the 
prompted  questions  provide  the  additional  information  needed  to  fill  the  remaining  data 
fields.  This  is  intended  to  be  a  first  step  towards  allowing  a  user  to  specify  or  design  an 
object  and  then  be  able  to  calculate  the  scattered  fields  from  this  object.  Although  it  is 
far  from  efficient,  it  does  serve  a  definite  purpose.  Further  refinement  of  this  program 
will  allow  for  an  easier  user  interface.  This  should  enable  technically  trained  personnel, 
without  a  detailed  understanding  of  these  programs,  to  benefit  from  the  Field  Feedback 
Formulation. 
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Figure  24.      Method  5  Mesh  Structure  Example 


Figure  25.      Method  6  Mesh  Structure  Example 
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Table  I.     SUMMARY  OF  MESH  GENERATION  CAPABILITIES 


Method 
Number 

Straight 

Line 
Bisection 

Straight 
Line  From 

Left  to 
Right  Side 

User  Se- 
lected Stop 
Node 

Equally 

Spaced 

Bisection 

Nodes 

User  Se- 
lected Start 
Node 

1 

X 

X 

■) 

X 

X 

X 

4 

X 

X 

5 

X 

X 

6 

X 

X 

X 
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IV.     FINITE  ELEMENT  BOUNDARY  VALUE  PROGRAM 

A.     INTRODUCTION 

The  Finite  Element  Boundan.^  Value  (FEBV)  program  is  the  feed  forw^ard  (U)  oper- 
ator in  the  Field  Feedback  Formulation  {F  )  as  shown  in  Figure  2  on  page  4.  This 
program  solves  the  Helmholtz  equation,  as  discussed  in  Chapter  II,  for  the  Dirichlet 
boundar\'  condition  specified  in  the  input  data  file,  as  discussed  in  Chapter  III.  These 
boundary  conditions  are  imposed  on  the  offset  boundary  contour.  The  purpose  of  this 
program  is  to  find  the  unknown  field  values  inside  the  offset  boundarv'  contour  within 
the  input  object.  Since  the  ultimate  goal  is  to  obtain  the  scattered  far  fields  from  the 
object,  the  unknown  field  values  on  the  perimeter  are  of  primar\'  interest.  It  will  also 
be  necessary  to  approximate  the  normal  derivative  of  the  field  at  the  object  perimeter. 
As  derived  in  Chapter  II.  the  goal  of  this  program  is  to  solve, 

The  A  matrix  represents  the  effect  that  the  I  —  1  row  has  on  the  /"'  row  values.  Similarly, 
the  B  matrix  represents  the  effect  that  the  I-th  row  has  on  itself,  while  the  C  matrix  re- 
presents the  effect  that  the  I  +  1  row  has  on  the  I-th  row.  All  other  rows  do  not  effect 
the  I-th  row.  This  is  due  to  the  pyramidal  basis  functions,  discussed  in  Chapter  II,  that 
all  have  zero  value  when  a  node  is  not  directly  connected  to  another  node.  The  P  matrix 
represents  the  combined  effects  of  the  boundary  nodes  on  the  I-th  row.  These  values 
are  transposed  to  the  right  side  of  the  equality  to  form  the  system  forcing  function.  It 
is  worth  remembering  that  for  the  I  =  1  row,  the  A  matrix  =  0  and  that  for  the  I  = 
I  MX  row,  the  C  matrix  =  0.  Thus,  using  the  row  by  row  stepping  process,  first  dis- 
cussed in  Chapter  III  of  the  mesh  generation  process,  the  A,  B,  C  and  P  matrices  can 
be  filled.  There  is  an  A,  B,  C  and  P  matrix  for  each  row,  and  it  is  necessar}^  to  calculate 
the  functional  for  two  elemental  rows  to  fill  one  set  of  matrices.  It  is,  therefore,  obvious 
that  the  data  is  used  twice,  once  for  the  current  row  and  again  for  the  next  row.  Spe- 
cifically, the  functionals  necessar\'  to  complete  the  B  matrix  and  totally  fill  the  C  matrix 
is  used  again  to  totally  fill  the  next  row's  A  matrix  and  partially  fill  the  B  matrix.  Thus, 
for  a  row  I  -  1,  (such  that  I  -  1  9^  I),  all  of  the  A  matrix  and  part  of  the  B  and  P  matrices 
are  filled.  When  the  row  is  stepped  to  I,  the  B  and  P  matrices  are  completed  and  all  of 
the  C  matrix  is  filled.    The  A,  B,  C  and  P  matrices  represent  tri-block  matrices  within  a 
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much  larger  global  matrix.  This  row  of  tri-block  matrices  in  the  global  matrix  can  be 
partially  solved  using  the  forward  portion  of  the  Ricatti  transform.  The  Ricatti  trans- 
form is  a  numerical  technique  that  optimizes  the  solution  for  banded  (tri-block)  systems 
of  linear  equations.   The  equations  are, 

and 

S,^,  =(5,  =  (5,  +  ^,/?,)-^(/>-4.S,) 

where  R,^^  is  the  i""  +  1  R  matrix  and  the  S,.,  is  the  /'''  +  1  S  vector.  Note  that  the  next 
rows  R  matrix  and  S  vector  depends  on  the  previous  rows  R  matrix  and  S  vector. 
During  the  for\vard  step  (MARCH  subroutine)  an  R  matrix  and  S  vector  is  generated 
for  each  row.  When  all  rows  have  been  calculated,  a  back  sweep  computes  the  unknown 
field  values,  »//,,.    The  equation  is 

lAz-i  =  Ri^i  +  S^. 

This  back  sweep  must  read  the  RS  data  from  the  disk  backwards.  This  is  a  storage  in- 
tensive process  for  all  but  the  most  trivial  problem.  The  field  values  iA,_,  is  first  found 
by  remembering  that  the  last  row  (  I  =  IMX  ).  C;v,a-  =  0.  With  C,vfA'  =  0  ,  R,_mx^\  =0 
and  ^,Mx  =  SiMx-\  '  which  is  the  last  calculated  S  vector.  The  recursion  continues  until 
all  of  the  ip,  's  have  been  calculated. 

B.     FINITE  ELEMENT  BOUNDARY  VALUE  PROGRAM 

The  eight  subroutines  comprising  the  FEBV  program  are  discussed  below. 

1.  Zero 

This  subroutine  fills  all  the  array  positions  of  the  A,  B,  C  and  P  matrices  with 
0  +  jO.  This  zeroing  is  necessary  after  all  calculations  concerning  a  given  row  are 
completed. 

2.  Varint  (Variational  Integration) 

This  subroutine  calculates  the  complex  functional  for  a  given  input  element. 
The  functionals  reflect  the  efiect  that  each  node  has  on  the  three  nodes  associated  with 
an  element.  The  subroutine  must  be  provided  with  the  X,  Y  coordinates  of  the  element 
nodes,  and  the  material  parameters  £,  and  ijl,.  These  functionals  are  returned  in  a  3  x  3 
complex  matrix.    The  area  of  the  input  element  is  also  provided  as  a  by-product  of  this 
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numerical  integration.  The  areas  of  each  element  are  sorted  and  used  to  find  the  largest 
and  smallest  elements.  Once  found,  the  largest  and  smallest  areas  are  combined  to  form 
the  area  ratio,  which  is  defined  as, 

maximum  area 
area  ratio  =  — — . 


minmium  area 


This  ratio  should  be  kept  as  small  as  possible.  For  a  circular  cylinder,  values  slightly  less 
than  2.0  are  possible.  For  more  complicated  objects,  the  ratio  can  get  considerably 
larger.  A  ratio  >  2.5  causes  a  screen  message  of  "YOU  SHOULD  CONSIDER 
ABORTING  THIS  RUN  AND  LOOKING  AT  THE  MESH  IN  CURVE  DIGITIZER. 
A  BETTER  METHOD  MAY  BE  AVAILABLE".  It  may  or  may  not  be  possible  to 
obtain  an  area  ratio  <  2.5.  The  intention  of  the  area  ratio  is  to  insure  that  a  uniform 
mesh  is  constructed  prior  to  attempting  a  problem  solution.  Smaller  area  ratios  are  in- 
dicative of  mesh  structures  that  do  not  have  grossly  difierent  element  sizes.  This  will 
lead  to  more  accurate  finite  element  solution. 

3.  Fill 

This  subroutine  calculates  and  stores  all  of  the  functional  for  an  entire  ele- 
mental row.  This  is  accomplished  by  repeated  VARINT  calls  for  each  element  that 
comprises  an  elemental  row. 

4.  BNDC  (Boundary  Condition) 

This  subroutine  calculates  and  stores  the  boundary  conditions  desired  for  the 
olTset  boundarv'  contour.  These  conditions  can  be  plane  wave,  cylindrical  modes  or  in- 
put manually.  For  plane  wave  conditions,  the  percent  error  of  the  FEBV  program  will 
also  be  returned.  This  calculation  only  has  meaning  for  the  uniform  case  (Field  5  set  to 
"U"  or  "u").  Memory  limitations  do  not  allow  this  feature  for  the  cylindrical  mode 
boundary  conditions  for  modes  other  than  zero  and  one.  A  separate  routine  could  be 
made  available  for  offline  percent  error  calculations.  No  such  feature  is  provided  if  the 
boundary  conditions  are  manually  input.  Any  desired  boundary  condition  may  be  gen- 
erated by  modifying  or  appending  the  proper  code  to  this  subroutine. 

5.  Loader 

This  subroutine  loads  the  A,  B,  C  and  P  matrices  for  a  given  row.  For  all  but 
the  I  =  1  and  IMX  rows,  two  calls  of  FILL/VARINT  for  two  elemental  rows  of  data 
are  required  to  fill  the  A,  B.  C  and  P  matrices.  This  is  an  additive  process  that  starts 
after  the  ZERO  subroutine  initializes  the  matrices.  Successive  functional  are  added  to 
the  existing  data  in  the  proper  array  positions. 
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6.  March 

This  subroutine  performs  the  forward  portion  of  the  Riccati  transform.  The 
output  data,  called  RS.DAT  (R  matrix  and  S  vector),  are  stored  on  disk. 

7.  CSMINV  (Complex  Square  Matrix  Inversion) 

This  subroutine  accomplishes  the  complex  matrix  inversion  required  by  the 
foru'ard  Riccati  transform.  A  maximum  dimension  of  50  x  50  was  established  to  limit 
memory  utilization.   This  allows  for  a  maximum  unknown  width  of  50. 

8.  Sweep 

All  of  the  above  subroutines  are  called  at  least  once  for  each  row.  The  sweep 
routine  is  called  only  after  all  of  the  row  calculations  are  completed,  and  then  only  once. 
This  subroutine  conducts  the  Riccati  back  sweep  by  reading  the  RS  data  generated  by 
the  MARCH  subroutine.  This  data  is  read  olT  the  disk  backwards,  using  the 
FORTRAN  "backspace"  command.  The  returned  field  values  are  actually  individual 
contributions  due  to  a  unit  valued  basis  function  being  individually  applied  to  each  of 
the  boundary  nodes.  In  so  doing,  the  problem  need  only  be  solved  once.  After  the  data 
is  stored,  in  matrix  form  as  the  U.DAT  file,  the  boundar>'  value  problem  may  be  solved 
for  any  incident  field  by  multiplying  the  U  matrix  by  the  new  incident  fields.    Thus. 

y^ perimeter  ~  L^'J  *  Y 'mcidenf 

A  summary  of  the  input,  output  and  error  data  is  provided  in  the  OUTPUT. DAT  file. 

9.  Save 

This  subroutine  was  necessarv"  to  allow  for  reasonably  sized  problems.  Since 
all  of  the  Field  Feedback  Formulation  code  could  not  fit  into  640  kilobytes  of  memor>', 
the  program  was  divided  into  two  parts.  All  of  the  data  necessary  to  perform  the  pro- 
grams is  saved  in  the  F3.DAT  file.  This  data  is  then  read  by  the  field  feedback  program. 
This  technique,  though  very  inefficient,  circumvents  the  640  kilobyte  memory  limitation 
imposed  by  the  IBM  Disk  Operating  System  (DOS).  EfTorts  to  convert  this  code  to  run 
under  a  compiler  that  does  not  have  a  640  kilobyte  memor\'  limitation,  such  as  Micro- 
way  NDP  FORTR^^N  compiler,  will  be  pursued  at  a  later  date. 

C.     VARINT  VALIDATION 

The  first  step  in  the  program  validation  required  an  understanding  of  the  error 
convergence  of  the  variational  elements  as  a  function  of  element  size,  d  and  material 
properties,  c,  and  u,  .  The  test  program  provided  in  Appendix  D  varies  the  element  size, 
in  wavelengths,  for  a  test  mesh  structure.    This  test  structure  shown  in  Figure  26  on 
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page  46  has  only  one  unknown  at  the  center.  The  mesh  is  made  of  four  adjacent  ele- 
ments. These  elements  are  mside  a  uniform  material  having  properties,  £,  and  /i,.  The 
other  four  nodes  are  established  as  boundary  nodes. 

PLANEWAVE 


Figure  26.      Test  Mesh  Structure 

A  plane  wave,  of  user  specified  frequency  and  amplitude,  is  used  to  determine  the 
boundar}'  conditions  on  the  four  boundary  nodes.  This  plane  wave  is  propagating  down 
the  vertically  oriented  axis.  The  unknown  nodal  value  is  calculated  and  compared  to  the 
actual  plane  wave  value  (1  +  jO).  A  percent  error  is  calculated  as  the  dimension  of  the 
elements,  d  are  reduced.  As  expected,  the  solution  became  more  accurate  as  the  ele- 
ments become  smaller.  A  plot  of  the  results  of  this  test  program  is  provided  in 
Figure  27  on  page  47.   This  data  was  generated  with  £,  =  1  +yO,  and  fj.,  =  1  +yO. 
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Figure  27.      Solution  Error  for  a  Test  Mesh  Structure 


The  region  of  interest  is  where  the  error  approaches  zero.  An  expanded  plot  of  this  re- 
gion is  provided  in  Figure  28  on  page  48.  1  he  accuracy  of  the  solutions  were  desired 
within  1  percent  of  the  exact  value.  This  requires  an  element  that  is  <"7y  .  in  the  ma- 
terial. To  ensure  the  desired  error  was  achieved,  elements  were  usually  scaled  to  be 
<  -rj-  ,  in  the  material.  The  phase  of  the  plane  wave  was  also  varied  to  determine  the 
eifect  on  convergence.   No  significant  effect  was  noted. 
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Figure  28.      Solution  Error  for  a  Test  Mesh  Structure  (Expanded) 

D.     FINITE  ELEMENT  BOUNDARY  VALUE  PROGRAM  VALIDATION 

1  he  next  validation  step  required  an  actual  object  to  calculate  the  fields  in  and  on. 
The  simplest  object  was  actually  not  an  object  at  all,  but  rather  an  imaginary  circular 
cylinder  in  free  space.  A  plane  wave  was  propagated  through  this  free  space.  Since  no 
material  interface  existed,  the  exact  solution  would  be  known  for  all  positions.  The 
plane  wave  established  the  boundary  conditions  on  the  offset  contour.  Trial  runs  were 
conducted  varying  the  mesh  resolution  and  boundary  offset  contour  distance.  With  the 
error  defined  as, 


error  = 


I(calculatcd  value  —  actual  value) 


I(actual  value)' 


two  errors  were  defined.  The  perimeter  error  only  considers  the  perimeter  nodes  in  the 
error  calculation.  The  bisection  segment  error  only  considers  the  bisection  segment 
nodes,  with  the  exception  of  the  two  end  nodes,  in  the  error  calculation.   The  end  nodes 
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of  the  bisection  segment  are  part  of  the  perimeter  and  are,  therefore,  not  considered. 
As  expected,  the  perimeter  error  could  be  rapidly  reduced  by  decreasing  the  offset  dis- 
tance. 1  his  is  due  to  the  fact  that  the  node  or  nodes  closest  to  an  unknown  node  dom- 
inate the  field  contributions  at  this  node.  Again,  a  goal  of  <  1%  error  was  desired.  1  he 
reduced  offset  distance  had  a  very  small  effect  on  the  bisection  segment  error.  The  only 
way  to  significantly  reduce  this  error  was  to  increase  the  mesh  resolution.  An  increased 
mesh  resolution  reduced  both  errors.  Figure  29  shows  the  perimeter  error  as  a  function 
of  the  number  of  bisection  segment  nodes. 


I 


Numbar  of  BItaetlon  Nod«« 


Figure  29.      Perimeter  Error 


The  three  curves  are  for  offset  distances  of  0.05,  0.025  and  0.01  X  .  Figure  30  on  page 
50  shows  the  bisection  segment  error  as  a  function  of  the  number  of  bisection  segment 
nodes.  The  two  curves  are  for  offset  distances  of  0.05  and  0.025  X^ .  For  offset  distances 
less  than  0.025  ).^  ,  the  curves  are  almost  identical  to  the  0.025  X^  curve.  These  curves 
are  omitted  for  clarity. 
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Figure  30.      Bisection  Segment  Error 


The  calculated  and  exact  field  values  for  a  0.5  Aq  diameter  circular  cylinder,  with  a 
0.05/^0  mesh  resolution  and  olTset  distance  and  £,  =  1  4-^0  is  shown  in  Figure  31  on  page 
51.  The  exact  fields  values  for  the  real  and  imaginary  portion  of  the  plane  wave  are 
shown  as  squares  and  diamonds,  respectively.  1  he  solid  curves  are  the  calculated  field 
values.  The  perimeter  and  bisection  errors  were  0.74  and  1.69  percent,  respectively. 
Figure  32  on  page  52  shows  the  elTects  of  not  properly  adjusting  the  mesh  resolution 
and  offset  distance  when  the  material  is  changed.  In  this  case,  the  permittivity  was 
changed  to  e,  =  4+J0.  Ihis  caused  the  perimeter  and  bisection  segment  errors  to  in- 
crease to  18.2  and  41.1  percent,  respectively.  The  mesh  resolution  and  ofiset  distance 
were  reduced  to  0.021/^.0  and  the  permittivity  was  changed  to  £,  =  4  —j4.  The  results  are 
shown  in  Figure  33  on  page  53  .  This  proper  selection  of  mesh  resolution  and  offset 
distance  reduced  the  perimeter  error  and  bisection  segment  error  to  0.44  and  0.56  per- 
cent, respectively.    Note  that  in  the  lossy  material  case  the  two  errors  are  very  close  in 
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magnitude.    This  is  due  to  the  way  that  the  error  is  defined.    This  definition,  in  a  lossy 
material,  overemphasizes  the  objects  leading  edge. 
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Figure  31.      0.5  ).  cylinder,  e,  =  1  +  jO 

This  portion  of  the  object  (as  well  as  the  trailing  edge)  is  the  most  accurate  for  the 
bisection  segment  calculations.  This  is  because  of  the  relative  closeness  of  these  nodes 
to  the  perimeter.  For  the  lossless  material,  the  bisection  segment  error  is  typically  two 
to  three  times  the  perimeter  error,  for  the  same  number  of  bisection  nodes.  The  con- 
clusion that  the  offset  distance  should  be  made  as  small  as  possible  in  order  to  minimize 
the  perimeter  error  is  correct,  but  will  prove  to  be  counterproductive  in  the  long  run. 
There  is  a  competing  elTect  that  will  require  this  contour  offset  distance  to  be  as  large 
as  possible.  This  effect,  associated  with  the  accuracy  of  the  Green's  function  contour 
integral,  will  be  discussed  in  Chapter  V. 
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E.     INHOMOGENEITY 

Until  this  point,  all  testing  was  done  in  circular  homogeneous  dielectric  materials. 
It  was,  at  a  minimum,  desirable  to  place  the  object  in  a  vacuum  to  calculate  the  scattered 
fields.  This  requires  the  first  two  and  last  two  elements  of  each  row  to  have 
c,  =  1  +/)  and  n,=  1  +jO.  These  elements  represent  the  space  between  the  objects  per- 
imeter and  the  ofTset  boundary  contour.  Since  the  plane  wave  solution  is  no  longer  valid 
for  this  geometry,  a  new  problem  with  a  known  solution  was  needed. 

Given  a  homogeneous  dielectric  circular  cylinder  of  radius  a  and  permittivity  c,  ,  as 
shown  in  Figure  34  on  page  54,  it  can  be  shown  that  [Ref  8  ], 

>A„(/^,  0)  =  A„J„{k,R)  cos  «0, 
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Figure  34.      Cylindrical  Mode  Geometry 


where, 


►        X 


ijy^iR,  (p)  is  the  field  value  inside  the  dielectric  material. 


"         Hl^\R,)lUkMJ',{R„)  -  krI'.ikMUR,)} 


Jn  =  Bessel  Function  of  the  l"  kind  of  order  n 


and 


K  =  ylh 
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An  additional  formulation  is  required  for  the  fields  outside  of  the  cylinder.  A  cylindrical 
wave  (mode,  n  =  1)  was  applied  to  a  homogeneous  dielectric  cylinder.  The  exact  and 
finite  element  field  values  were  calculated  and  compared.  Figure  35  on  page  55  ,  shows 
a  typical  result.  Mesh  resolution,  olfset  distance  and  material  permittivity  were  varied 
to  determine  convergence,  1  he  results  were  similar  to  the  dielectric  homogeneous  plane 
wave  case  discussed  earlier. 


ParimaUr  Noda  Numbar 

Figure  35.      Typical  Cylindrical  Mode  Boundary  Value  Problem  Result 


F.     FINITE  ELEMENT  CONCLUSIONS 


High  accuracy  solutions  can  be  obtained  by  maintaining  a  mesh  resolution  of 
<-r-r  in  the  material.    The  offset  distance  should  be  kept  at  approximately  the  same 


;. 


magnitude  as  the  mesh  resolution,  but  mav  be  as  large  as  -r— .  Increased  accuracy  re- 
quires  an  increased  mesh  resolution.  This  increase  in  accuracy  is  at  the  expense  of  in- 
creased processor  time. 


55 


V.     FIELD  FEEDBACK  PROGRAM 

A.     INTRODUCTION 

The  Field  Feedback  Program  is  the  feedback  (T)  operator  in  the  F .  This  program 
utilizes  the  output  of  the  finite  element  program  and  the  input  incident  fields  to  evaluate 
a  "near  field"  Green's  function  integral.  This  integral  is  called  "near  field"  in  that  the 
integral  will  be  performed  within  -r—  of  the  object  perimeter.    As  seen  in  Figure  36, 


three  surface  contours  are  defined, 
boundary-  condition  is  applied. 


The  boundarv  contour  is  where  the  incident  field 


Geometric  Contour,  GC 


Boundary,  B 


Figure  36.      Contour  Arrangement 


The  perimeter  contour  is  where  the  finite  element  program  solves  for  the  field  values  on 

the  object  perimeter.    The  geometric  contour  is  midway  between  the  boundary  and  the 

perimeter  and  is  the  contour  over  which  the  Green's  function  contour  integral  is  per- 

c  ^u 
formed.    This  is  necessary  to  allow  the  — —  to  be  approximated  by  a  finite  difibrence 
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technique.    The  \l/cc  is  the  average  of  the  field  values  on  the  boundary  and  perimeter. 
Thus, 

dp  yy  boundary  ~  r  perimeter) 

and 

y^f^ perimeter   '    r boundary' 


^GC  = 


The  Field  Feedback  program  is  provided  as  Appendix  E. 

B.     FIELD  FEEDBACK  PROGRAM 

1.  Input 

This  subroutine  reads  the  finite  element  data  stored  in  the  F3.DAT  file  by  the 
SAVE  subroutine.   All  necessary  nodal  X.  Y  coordinates  are  calculated. 

2.  TMAT  (T  Matrix) 

This  subroutine  loads  a  complex  matrix  called  TMAT.  Each  element  of  the  T 
matrix  is  the  result  of  a  Green's  function  contour  integral.  This  integral  is  conducted 
on  the  GC  contour.  The  m-th  column  of  the  T  matrix  is  evaluated  with  a  single  unit 
valued  basis  function  on  the  m-th  boundary  node.   The  equation, 


T       = 

"■  n,  m 


^^n  .       CG 

-^ ^n^— 

cn  en 


'GC 


dGC, 


may  now  be  numerically  evaluated  at  each  of  the  n  boundar}'  nodes.    This  process  is 
repeated  for  each  row  until  the  entire  matrix  is  filled. 

The  numeric  integration  uses  a  rectangular  mid-point  approximation  technique. 
For  this  linearized  problem,  this  is  equivalent  to  a  trapezoidal  integral  technique.  The 
integral  path  between  any  given  two  nodes  is  subdivided  based  on  the  distance  to  the 
first  point.  This  subdivision  maybe  controlled  by  three  input  variable  fields.  A  multi- 
plicative scale  factor  and  offset  terms  are  available.  To  date,  the  utilization  of  the  scale 
factor  (set  >  1.0)  and  the  offset  terms  (set  >  0  )  have  not  proved  necessar}'.  This  may 
be  necessar\"  for  more  complicated  geometric  structures. 
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3.     CNSOLV 

This  subroutine  solves  the  equation, 


Cn  =  U  -  T„,amx'\       '  ^ 


boundary^ 


where  C„  is  the  total  scattered  fields  on  the  boundary',  I  is  the  identity  matrix  and 
^boundcry  is  the  initially  specified  incident  fields. 
4.     FFLD  (Far  Fields) 

This  subroutine  calculates  the  scattered  far  field  or  bistatic  radar  cross  section 
per  unit  length  of  the  object.  Any  desired  integer  angular  resolution  for  the  calculation 
may  be  specified.   The  final  results  are  stored  in  the  FFPAT.DAT  file. 

C.     FIELD  FEEDBACK  VALIDATION 

Given  a  plane  wave  boundary  conditions  imposed  on  an  imaginar\'  cylinder  in  free 
space,  a  Green's  function  contour  integral  may  be  performed  around  the  cylinder.  This 
cyhnder  is  imaginary  in  the  sense  that  it  does  not  actually  exist.  The  cylinder  is  actually 
an  artificial  boundar\'  that  does  not  form  a  material  interface.  For  points  inside  this 
cylinder,  the  resulting  integral  should  equal  the  negative  of  the  actual  plane  wave  value 
at  that  point  in  freespace.  For  points  outside  the  cylinder,  the  resulting  integral  should 
equal  zero.  Several  test  conditions  were  investigated  to  ensure  that  the  numerical  inte- 
gration of  the  Green's  function  did  arrive  at  the  expected  values.    Maximum  absolute 

errors  for  these  tests  were  less  than  7  x  10"'  .    These  test  cases  started  with  exact  values 

ci// 
of  the  field.  \p  and  the  normal  derivative,  ^:^ —  .    After  initial  validation,  finite  element 

calculated  ip  and  — —  values  were  used.    Errors  increased  bv  approximatelv  a  factor  of 
en 

10.    Numerous  competing  effects  were  observed  with  two  dominant  effects  becoming 
evident. 

1.     Small  Object  Phenomenon 

For  very  small  objects,  where  only  a  few  elements  are  needed  to  meet  the  -rr- 
mesh  resolution  requirement,  the  normal  derivative  accuracy  is  crucial.  Since  the  normal 
derivative  is  the  calculated  normal  to  the  piecewise  linear  approximation  of  the  objects 
perimeter,  this  fit  is  usually  not  adequate.  To  prevent  this  problem,  additional 
perimeter/boundary  nodes  are  needed.  This  is  easily  accomplished  by  increasing  the 
mesh  density,  which  is  controlled  by  the  mesh  resolution  (input  field  7).  For  these  ver\' 
small  objects,  the  RCS  is  almost  uniform  for  the  T.M  case  and  co-sinusoidal  for  the  TE 
case.   The  utility  of  this  calculation  must,  therefore,  be  questioned. 
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2.     Offset  Distance  Phenomenon 

The  ofTset  distance  must  not  be  made  too  small  or  the  accuracy  of  the  Green's 
function  integral  will  sufTer.  For  a  circular  cylinder,  the  optimal  offset  was  between  0.03 
/  and  0.035  /.   Note  that  these  distances  are  alwavs  < 


/. 


20 
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VI.     VALIDATION 

A.     INTRODUCTION 

The  difiiculty  in  the  total  Field  Feedback  Formulation  {F  )  program  validation  is  in 
finding  problems  that  have  well  established  or  accepted  solutions.  The  total  F  program 
is  the  combination  of  the  mesh  generation.' finite  element  program  and  the  field  feedback 
program.    If  possible,  a  problem  that  has  a  closed  form  analytical  solution  is  desired. 


B.     HOMOGENEOUS  CIRCULAR  CYLINDRICAL  SCATTERING 

This  is  the  first  test  case  for  the  F  program.  A  penetrable  circular  cylinder  has  a 
closed  form  analytic  solution,  so  the  final  results  could  be  verified  against  exact  values. 
It  can  be  shown  that  for  the  TM  case  (E  -  wave)  the  reflection  coefficient  is, 


^  n 


UkMJ',{R,)  -  J-^  J',{k,R,)J„{Ra) 


UkMH^^^'iR,)  -  J-^  J',{k,RM\Ra) 


and  that  for  the  TE  case  (H  -  wave)  the  refection  coefficient  is, 


pTM  _ 


Uk,R,)J'„{R,)  -     1^  J'.ikMURa) 


V    «r 


r(2)'/ 


Uk,R,)Hl;>  {RJ  -  J^  J',ikMHl^\R,) 


,/Z 


where  k,  =  J h.e,  ,     y\  =  yj  -j^       and      2,=  ^ -^  .  [Ref  9  ]  The  scattered  field  for  the 
TM  case  is. 


y". 


4(7?,  0)  =  -  £' 


incident 


r^H'i\R)  +  2^;-T„//?\/?)  cos  n4> 


n=\ 


Similarly,  the  scattered  field  for  the  TE  case  is, 


//J(/?,  4>)  =  -H 


incident 


00 

V,H^^\R)  +  2^y- "r„//f  (;?)  cos  n<i> 


n=\ 
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In  the  far  field  region,  as  R  approaches  infinity,  the  scattering  width  may  be  defined  as. 


o{(t))  =  In  lim  p 


W\' 


p-H.00  •       I  77incident  1  2 


and 


(^)  =  y-lro  +  2^r,cosA70r 


The  source  code,  without  the  Bessel  fiinction  routines,  for  these  calculations  is  provided 
as  Appendix  F. 

Three  different  radius  dielectric  cylinders  were  tested.  Wavenumber  normalized  radii 
of  0.5,  1.0  and  2.0  with  c,  =  2.56  +  jO  were  selected.  The  TE  and  TM  results  of  the  0.5, 
1.0  and  2.0  radii  problems  are  provided  as  Figure  37  on  page  62,  Figure  3S  on  page 
63,  and  Figure  39  on  page  64  .  To  validate  the  Chapter  VI  conclusion  that  the  unit 
normal  was  the  cause  of  the  errors  for  ver\'  small  circular  cylinders,  the  £,  was  increased 
from  2.56  to  25.6.  This  value  still  meets  the  requirement  for  mesh  resolution  not  to  ex- 
ceed  -r^  .  The  TE  and  TM  results  are  provided  as  Figure  40  on  page  65  .  Actual  cross 
section  values  are  indicated  by  the  squares  and  the  Field  Feedback  Formulation  sol- 
utions are  plotted  as  a  single  solid  curve. 

C.     HOMOGENEOUS  IRREGULAR  OBJECTS 

The  results  detailed  in  [Ref  10]  were  next  validated  with  the  F  program.  This  ob- 
ject, as  seen  in  Figure  41  on  page  66  ,  is  a  dielectric  shell  with  inner  radius  =  .25  /g, 
outer  radius  =  .30  /q  and  £,  =  4  +  jO.    The  comparison  of  the  [Ref  10]  data  and  the 

3 

f  results  are  provided  as  Figure  42  on  page  67. 
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PLANEWAVE 


Figure  41. 


Dielectric  Shell  Mesh 


The  TE  case  shows  excellent  agreement.    The  TM  case  tends  to  diverge  at  the  smaller 
values  of  0. 

D.     AN  INHOMOGENEOUS  OBJECT 

The  results  detailed  in  [Ref.  11]  were  next  validated  with  the  F  program.  This  ob- 
ject, as  seen  in  Figure  43  on  page  68,  is  a  dielectric  ring  with  inner  radius  =  .25  X^ ,  outer 
radius  =  .30  X^  and  e,  =  4  +  jO.  7  he  exact  solution  to  this  problem  is  available.  [Ref. 
11] 
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Figure  43.      Dielectric  Ring 

The  mesh  generation  program  was  not  modified  to  conform  to  both  the  outer  and  inner 
material  boundaries.  The  original  mesh  generation  program  design  concept  was  to 
closely  fit  the  mesh  to  the  objects  outer  perimeter  and  then  allow  for  material  inhomo- 
geneities  to  be  accounted  for  by  difierent  material  parameters  being  assigned  to  individ- 
ual elements.  A  section  of  the  generated  mesh  is  shown  in  Figure  44  on  page  69.  Note 
the  additional  curved  line  showing  the  .25  }.^  inner  radius.  This  curve  does  not  follow 
the  established  element  boundaries.  The  efl'ective  ring  is  shown  in  Figure  45  on  page 
69.  This  resulted  in  the  inner  radius  having  an  irregular  pattern  as  elements  near  the 
material  interface  varied  in  material  composition.  This  is  similar  to  the  granular  noise 
problem  characteristic  of  delta  modulation  communication  channels.  The  TE  and  TM 
results  are  provided  as  Figure  46  on  page  70. 
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Figure  44.      Partial  Mesh  nitli  Inner  Radius  Curve 


Figure  45.      Effective  Geometry  for  the  Dielectric  Mesh 
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VII.     CONCLUSIONS 

A.  RESULTS 

The  Field  Feedback  Formulation  proved  to  be  an  excellent  tool  for  calculating  the 
internal  and  scattered  fields  fi-om  the  tested  inhomogeneous  asymmetric  objects.  The 
keys  to  satisfactory  results  are, 

•  Keep  the  maximum  element  dimension  <  -rr-  , 

•  Maintain  a  mesh  with  a  uniform  distribution,  and 

•  Maintain  the  offset  distance  near  the  mesh  resolution  in  magnitude,  but  no  greater 
than  ^  . 

The  techniques  used  proved  to  be  capable  of  adapting  to  a  wide  variety  of  situations. 
These  adaptations  did  not  require  program  modifications  or  reprogramming.  The  nu- 
merous built-in  features,  as  discussed  in  the  input  field  description  of  Chapter  III,  al- 
lowed for  this  robustness. 

B.  RECOMMENDATIONS  AND  EXTENSIONS 

With  the  basic  program  validated,  future  testing  and  development  should, 

•  Emphasize  large  object  vahdation  against  accepted  results. 

•  Stress  inhomogeneity  and  irregularity  in  all  testing. 

•  Optimize  the  T  matrix  element  calculation  by  improving  the  Green's  function 
contour  integral. 

•  Evaluate  the  usefulness  of  the  maximum  distance  feature  (Field  12).  Beyond  this 
maximum  distance  no  contribution  is  made  to  the  Green's  function  integral.  It  is 
not  expected  that  this  will  be  possible  for  objects  less  than  a  few  wavelengths  in 
maximum  dimension. 

•  Modify  the  mesh  generation  program  to  allow  for  conformity  to  multiple  interfaces 
(multi-layered  objects  without  the  granular  noise  problem). 

•  Modify  all  programs  for  Intel  80386  based  Fortran  compiler  use,  such  as  NDP 
FORTRAN  by  Vlicroway.  This  will  remove  the  640  KB  memor\'  limitation  im- 
posed by  the  IBM  DOS  and  allows  for  the  solution  to  larger  scattering  problems. 
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APPENDIX  A.     INPUT  DATA  FILE  EXAMPLE 


CIRCLE 
R 
NI 
ND 
NU 
NM 

0.019 
0.03 
1.0 
0 
0 

999.  9 
36 
72 
1 
5 

35 
4.5 
4.0 
0.  3 
1.0 

0.  0 

1.  0 
0.  0 

p 

1.0 

3.  OOE+08 
0 

10 

20 

30 

40 

50 

60 

70 

80 

90 
100 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 


PLOT  TITLE 

RECTANGULAR  INPUT  DATA 

NO  INTERMEDIATE  DATA  RECORDING 

NO  DISSPLA  PROGRAM  GENERATION 

NON  UNIFORM  MATERIAL  (MATERIAL  INTERFACE  PRESENT) 

NO  MESH  GENERATION  ONLY 

REQUESTED  MESH  RESOLUTION  (IN  WAVELENGTHS) 

CONTOUR  DISTANCE  (IN  WAVELENGTHS) 

MULTIPLICATIVE  GFI  SCALE  FACTOR 

DISTANCE  <  1.  0  BIAS  TERM 

DISTANCE  >  1.  0  BIAS  TERM 

MAX  DIST  BEYOND  WHICH  THERE  IS  NO  CONTRIBUTION  TO  GFI 

NUMBER  OF  INPUT  DATA  POINTS 

NUMBER  OF  POINTS  FOR  SIGMA  CALCULATION  (IN  THE  CIRCLE) 

ELEMENT  GENERATION  METHOD 

START  NODE 

STOP  NODE 

X  AXIS  OFFSET 

Y  AXIS  OFFSET 

DESIRED  DIMENSION  (ORIGIN  TO  FIRST  POINT  .WAVELENGTHS) 

REAL  PART  OF  ALPHA 

PART  OF  ALPHA 

OF  BETA 

PART  OF  BETA 

GENERATOR  ENABLED 

AMPLITUDE 

PLANEWAVE,  FO 

-  ANGLE,  X,  Y  COORD. 


IMAGINARY 
REAL  PART 
IMAGINARY 
PLANEWAVE 
PLANEWAVE 
INPUT  FREQUENCY  OF  THE 


00000000 
08682409 
17101010 
25000000 
32139380 
38302220 
43301270 
46984630 
49240390 
50000000 
49240390 
46984630 
43301270 
38302220 
32139380 
25000000 
17101000 
08682404 
00000004 
08682413 
17101010 
25000000 
32139380 
38302220 


50000000 
49240390 
46984630 
43301270 
38302220 
32139380 
25000000 
17101010 
08682407 
00000002 
08682411 
17101010 
25000000 
32139380 
38302220 
43301270 
46984630 
49240390 
50000000 
49240390 
46984630 
43301270 
38302220 
32139380 
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240 

-.43301270 

-.  25000000 

250 

-.46984630 

-.  17101000 

260 

-.49240390 

-.08682403 

270 

-.50000000 

. 00000007 

280 

-.49240390 

. 08682416 

290 

-.46984630 

. 17101010 

300 

-.43301270 

. 25000010 

310 

-.  38302220 

.32139390 

320 

-.32139380 

. 38302230 

330 

-. 24999990 

.43301280 

340 

-. 17101000 

.46984630 

350 

-.  08682401 

.49240390 
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APPENDIX  B.  READ  PROGRAM 

REAL  A(0: 2000),  B(0:2000),  MRES,  DIST,  XORIGIN,  YORIGIN 

REAL  C,  D,  E,  F,  DPER,  DD,  FREQ,  EO,  MAXD 

INTEGER  I,  J,  METHOD,  STARTND,  STOPND,  NRES ,  MODE,  LBIAS,  GBIAS 

CHARACTER^^  1  CHAR ,  CHAR  1 ,  CHAR2 ,  CHARS ,  CHAR4 ,  CHARS 

CHARACTER^'-IS  NAME 

OPEN(UNIT  =  1,  FILE  =  'D:  FINALDWG.DAT') 
OPEN(UNIT  =  2,  FILE  =  'C:  MSFORT  INPUT.DAT') 
J  =  0 

VRITE(",'0  'ENTER  THE  PLOT  NAME  OR  LABEL  (MAX  OF  13  CHARACTERS)' 

READ(''-,1005)  NAME 

WRITEC^S^O  'ENTER  THE  COORDINATE  SYSTEM  IN  USE.  (R  OR  P)' 

READ(^'SlOOO)  CHAR 

WRITE(''S^--)  'DO  YOU  WANT  INTERMEDIATE  VALUES  ?  (I)' 

READ(^'--,1000)  CHARl 

WRITE(^S-'0  'DO  YOU  WANT  A  DISSPLA  PROGRAM  GENERATED  ?  (D)' 

READ(^'SlOOO)  CHAR2 

WRITE(--'-,''0  'UNIFORM  SLAB  (NO  MATERIAL  TO  VACUUM  INTERFACE  ?  (U)' 

READ(->,1000)  CHARS 

VRITE(^S-''-)  'MESH  GENERATION  ONLY  ?  (M)' 

READ(--'slOOO)  CHAR4 

WRITE(''-,^'0  'ENTER  THE  MESH  RESOLUTION.  (0.4,  ETC...)' 

READ(''S^'-)  MRES 

WRITE(^'S^'0  'ENTER  THE  CONTOUR  DISTANCE.  (O.S,  ETC...)' 

READ(-S-'0  DIST 

WRITE('S--'^)  'ENTER  THE  GFI  SCALE  FACTOR.  (1.0,  ETC...)' 

READ(-S^--)  DPER 

WRITE(''-,''0  'ENTER  THE  GFI  (  <  1  BIAS  TERM  (0,1,  ETC...)' 

READ( -'■■-,''-)  LBIAS 

WRITE (■-'%-'■-)  'ENTER  THE  GFI  (  >  1  BIAS  TERM  (0,1,  ETC...)' 

READ(-'S-'-)  GBIAS 

WRITE('-,-'0  'ENTER  THE  GFI  MAX  DISTANCE  (999.0,  ETC...)' 

READ(^>-,^0  MAXD 

WRITE('V,Vr)  'ENTER  THE  NUMBER  OF  POINTS  FOR  SIGMA  CALC.  ( S6 ,ETC.  .  .  ) ' 

READ('S^'-)  NRES 

WRITE(*,''0  'ENTER  THE  DRAWING  METHOD.  (1  -  6)' 

READ(^S^O  METHOD 

WRITE(*,''0  'ENTER  THE  STARTING  NODE  NUMBER.  (MUST  BE  >  0)' 

READ(--'-,^0  STARTND 

WRITE(*,''0  'ENTER  THE  BISECTION  STOPPING  NODE  NUMBER.  ' 

READ(-'S^O  STOPND 

WRITE('^''^)  'DESIRED  DISTANCE  FROM  ORIGIN  TO  POINT  1  (IN  WLs)  ' 

READ(^S^--)  DD 

WRITE('>,''0  'ENTER  THE  X  AXIS  ORIGIN.  ' 

READ(''-,^0  XORIGIN 

WRITE(->,''0  'ENTER  THE  Y  AXIS  ORIGIN.  ' 

READ(",^0  YORIGIN 

WRITE(-'S''0  'ENTER  THE  REAL  COMPONENT  OF  ALPHA.  ' 

READ('S-0  C 
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ENTER  THE  IMAGINARY  COMPONENT  OF  ALPHA. 


ENTER  THE  REAL  COMPONENT  OF  BETA. 


ENTER  THE  IMAGINARY  COMPONENT  OF  BETA. 


WRITEC^'S'"^) 

READ(''s-'0  D 

WRITE(''--,'0 

READ  (-'%■'-)  E 

WRITE(''s-''0 

RE  AD  (--■%■•'-)  F 

WRITEC^'S'O  'PLANE  WAVE  (P) 

READ(-'slOOO)  CHAR5 

WRITE C*,'"^)  'ENTER  THE  WAVE 

READ(''s-'0  FREQ 

WRITE( '•%'>)  'ENTER  THE  WAVE  AMPLITUDE 

READ(-'^*)  EO 

IF((CHAR5.EQ.  'C').OR.  (CHARS.  EQ 

WRITE ('>,*)  'ENTER  THE  MODE  NUMBER 

READ(*,''0  MODE 
END  IF 
DO  10,  I  =  0,  1999 

READdj^-sERR  =  20)   A(J),    B(J) 

IF(A(J).GT.  1000)  THEN 
GOTO  5 

ELSEIF(B(J).GT.  1000)  THEN 
GOTO  5 

ELSE 

J  =  J  +  1 


OR  CYLINDRICAL  MODE  (C) 
FREQUENCY  (IN  HERTZ)' 


c'))  THEN 


END  IF 

5     CONTINUE 

10    CONTINUE 

20    J  =  J  -  1 

WRITE(2,1005: 

)  NAME 

WRITE(2,1000: 

)  CHAR 

WRITE(2,1000: 

)  CHARl 

WRITE(2,1000: 

)  CHAR2 

WRITE(2,1000; 

)  CHARS 

WRITE(2,1000: 

)  CHAR4 

WRITE(2,1010: 

)  MRES 

WRITE(2,1010: 

)  DIST 

WRITE(2,1010; 

)  DPER 

WRITE(2,1020: 

)  LBIAS 

WRITE(2,1020~ 

)  GBIAS 

WRITE(2,1040; 

)  MAXD 

WRITE(2,1020; 

)  J 

WRITE(2,1020; 

)  NRES 

WRITE(2,1020' 

)  METHOD 

WRITE(2,1020' 

)  STARTND 

WRITE(2,1020; 

)  STOPND 

WRITE(2,1010: 

)  XORIGIN 

WRITE(2,1010' 

)  YORIGIN 

WRITE(2,1010' 

)  DD 

WRITE(2,1010' 

)  c 

WRITE(2,1010- 

)  D 

WRITE(2,1010 

)  E 

WRITE(2,1010 

)  F 

WRITE(2,1000 

)  CHARS 

WRITE(2,1010 

)  EO 

IF( ( CHARS. EQ. 

'C').OR.  (CHARS 

EQ. 


'c')) 


THEN 


WRITE (2, 1020)  MODE 
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END  IF 

WRITE ( 2, ''0  FREQ 

DO  30,  I  =  1,  J 

WRITE(2,1030)  I,  A(I),  B(I) 
30    CONTINUE 

WRITE(2,''0  'END  OF  FILE' 
C 

1000  FORMAT(Al) 
1005   F0RMAT(A13) 
1010  F0RMAT(F8.6) 
1020  F0RMAT(I3) 

1030     F0RMAT(1X,I3,7X,F12.8,5X,F12. 8) 
1040     F0RiMAT(F8.  3) 
C 

CLOSE(l) 

CL0SE(2) 

STOP 

END 
C 
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APPENDIX  C.     MESH  GENERATION/FINITE  ELEMENT  PROGRAM 

c 

C  FINITE  ELEMENT  MESH  GENERATION  PROGRAM 

C  WRITTEN  BY  LT.  T.  B.  WELCH 

C 

C  w/  PROGRAMMING  IDEAS  FROM  PROF.  M. A.  MORGAN 

C 

C  COMPLETELY  FILE  (INPUT.DAT)  DRIVEN 

C 

C  INPUT  PARAMETERS: 

C 

C  CHARACTER  (POLAR  OR  RECTANGULAR  INPUT  DATA  -  FLAG) 

C  CHARACTER  (INTERMEDIATE  MATRICES  WRITE  -  FLAG) 

C  CHARACTER  (DISSPLA  PROGRAM  GENERATION  -  FLAG) 

C  CHARACTER  (UNIFORM  MATERIAL  -  FLAG) 

C  CHARACTER  (MESH  GENERATION  ONLY  -  FLAG) 

C  MESH  RESOLUTION 

C  CONTOUR  OFFSET 

C  DESIRED  SCALE  FACTOR  FOR  GREEN'S  FUNCTION  INTEGRAL 

C  BIAS  (SHIFT)  FOR  GFI  STEP  WHEN  <  1. 0 

C  BIAS  (SHIFT)  FOR  GFI  STEP  WHEN  >  1.  0 

C  MAXIMUM  DISTANCE  BEYOND  WHICH  THERE  WILL  BE  NO 

C  CONTRIBUTION  TO  THE  GREEN'S  FUNCTION  INTEGRAL 

C  NUMBER  OF  POINTS 

C  NUMBER  OF  POINTS  FOR  CROSS  SECTION,  (IN  THE  CIRCLE) 

C  MESH  GENERATION  TECHNIQUE 

C  START  NODE 

C  STOP  NODE 

C  X  ORIGIN 

C  Y  ORIGIN 

C  DESIRED  DISTANCE  (FROM  ORIGIN  TO  FIRST  POINT,  WAVELENGTH) 

C  ALPHA  (INPUT  AS  A  AND  B  THEN  CONVERTED  TO  COMPLEX  ALPHA) 

C  BETA   (INPUT  AS  A  AND  B  THEN  CONVERTED  TO  COMPLEX  BETA) 

C  CHARACTER  (PLANEWAVE  GENERATOR  -  FLAG) 

C  AMPLITUDE 

C  FREQUENCY  (IN  HERTZ) 

C  --  OR  -- 

C  AMPLITUDE 

C  CYLINDRICAL  MODE  NUMBER  (GENERATOR  ALWAYS  DOES  N=l) 

C  FREQUENCY  (IN  HERTZ) 

C  --  OR  -- 

C  BOUNDARY  CONDITIONS 

C  DATA  POINT  PAIRS  (X,  Y  OR  RADIUS,  THETA) 

C  END  OF  FILE  MARKER 

C 

C 

COMPLEX  ALPHA, BETA, BCOND( 100), LINE(50),ANS( 100) 

COMPLEX  A(50,50),B(50,50),C(50,50),P(50,100),SURBC(100) 

COMPLEX  U(100,100),PSI(50,100),SVEC(50,100) 

REAL  MRES,MESH(0: 200,5) ,NDP(200 ,2) , OFFSET, DPER,MRESW 

REAL  MINAREA,MAXAREA,AREA,E0,XORIGIN,YORIGIN,K0,MAXD 

INTEGER  NPNTS,PERND,NODES(200),MOR(200),BIND,NBMX,UNK,LBIAS,GBIAS 
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INTEGER  LND( 0:  200 , 3) ,NDL( 200,4) ,NCT( 200) ,MAXEL, IMX ,MODE ,NRES 
INTEGER  METHOD , STARTND , STOPND , MINROW , MINEL , MAXROW , NABC( 100,3) 
CHARACTER--'--!  CHAR,  CHARl,  CHAR2 ,  CHAR3 ,  CHAR4,  CHAR5 
EQUIVALENCE  (PSI,  P) 

COMMON/BLKl/MESH 

C0MM0N/BLK2/PERND , BIND 

C0MM0N/BLK3/N0DES 

C0MM0N/BLK4/LND , NDL , NCT 

C0MM0N/BLK5/NDP 

C0MM0N/BLK6/M0R 

C0MM0N/BLK7/MINAREA, MINROW, MINEL, MAXAREA, MAXROW, MAXEL, AREA 

C0MM0N/BLK8/A,B,C,P 

C0MM0N/BLK9/CHAR2,  CHAR3,  CHAR4 


OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 

CACCESS 
OPEN  ( 

CACCESS 
OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 
OPEN  ( 


UNIT  =  1,  FILE  = 
UNIT  =  2,  FILE  = 
UNIT  =  3,  FILE  = 
UNIT  =  4,  FILE  = 
UNIT  =  7,  FILE  = 
=' SEQUENTIAL' ,FORM 
UNIT  =  8,  FILE  = 
=' SEQUENTIAL' ,FORM 
UNIT  =   9,  FILE  = 


UNIT  = 

10, 

FILE  =  ' 

UNIT  = 

11. 

FILE  =  • 

UNIT  = 

12, 

FILE  =  ' 

UNIT  = 

13, 

FILE  =  ' 

UNIT  = 

19, 

FILE  =  • 

UNIT  = 

20, 

FILE  =  ' 

UNIT  = 

30, 

FILE  =  ' 

UNIT  =  40,  FILE  = 


MATRIX.  DAT' , STATUS  =  'UNKNOWN') 
PLT. DAT' , STATUS  =  'UNKNOWN') 
INPUT. DAT' , STATUS  =  'OLD') 
TEXT. LBL' , STATUS  =  'UNKNOWN') 
RS. DAT' , STATUS  =  'UNKNOWN', 
'UNFORMATTED') 

PSI. DAT' , STATUS  =  'UNKNOWN' , 
'FORMATTED') 

ABCP. DAT' .STATUS  =  'UNKNOWN') 
BCOND.DAT',  STATUS  =  'UNKNOWN') 
FMATRIX. DAT' ,STATUS=' UNKNOWN' ) 
NABC. DAT' , STATUS=' UNKNOWN' ) 
OUTPUT. DAT ' , STATUS= ' UNKNOWN ' ) 
ABRMAT. DAT ' , STATUS= ' UNKNOWN ' ) 
DISPLA. DAT' ,STATUS=' UNKNOWN' ) 
U. DAT ' , STATUS= ' UNKNOWN ' ) 
F3. DAT' , STATUS= ' UNKNOWN ' ) 


MINAREA  =  999999. 9 
MAXAREA  =  -999999. 9 

CALL  I0( NPNTS , MRES , METHOD , STARTND , STOPND , OFFSET , ALPHA , BETA , BCOND , 
CEO, CHARl, MODE, XORIGIN,YORIGIN, CHARS, DPER,KO,NRES,MRESW,LBIAS,GBIAS 
C,MAXD) 

CALL  ROTATE ( NPNTS , METHOD , STARTND , STOPND ) 

CALL  B OUND ( MRE  S , NPNTS , METHOD , STOPND ) 

CALL  NORMAL 

CALL  NODSETC METHOD, NABC, NBMX,UNK) 

PAUSE  'PLEASE  PRESS  ENTER  TO  CONTINUE,  OR  CTRL  BREAK  TO  ABORT!  ' 

CALL  LOADERC  BCOND , OFFSET , ALPHA , BETA , NABC , IMX , NBMX , EO , SURBC , CHARl , 
CLINE, MODE, XORIGIN,YORIGIN,SVEC, CHARS) 

CALL  SWEEP( IMX , NABC , SURBC , LINE , CHARl , U , BCOND , PS I , ANS , CHARS ) 

CALL  SAVE( BCOND , ANS , U , OFFSET , PERND , CHARS , DPER , KO , XORIGIN , YORIGIN , 
CNRES,MRESW,LBIAS,GBIAS,MAXD) 

CLOSE(l) 
CL0SE(2) 
CL0SE(3) 
CL0SE(4) 
CL0SE(7) 
CL0SE(8) 
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CL0SE(9) 

CLOSE(IO) 

CLOSE(ll) 

CL0SE(12) 

CL0SE(13) 

CL0SE(19) 

CL0SE(20) 

CLOSE(30) 

CLOSE(40) 

STOP 
END 


SUBROUTINE  IO(NPNTS ,MRES , METHOD , STARTND , STOPND ,DI ST, ALPHA , BETA, 
CBCOND,EO,CHAR1, MODE, XORIGIN,YORIGIN, CHARS, DPER,KO,NRES,MRESW, 
CLBIAS,GBIAS,MAXD) 

THIS  SUBROUTINE  READS  IN  THE  INPUT  PARAMETERS 
AND  SURFACE  DATA  POINTS.  THESE  POINTS  CAN  BE 
IN  EITHER  POLAR  OR  RECTANGULAR  FORM. 

COMPLEX  ALPHA , BETA , BCOND( 100 ) 

REAL  A,B,PI,XORIGIN,YORIGIN,E0,K0,F0,SFAC,DD 

REAL  MESH(0: 200 ,5) ,MRES ,THETA, RADIUS ,DIST,C,LAMBDA,MRESW 

REAL  DISTW,DPER,MAXD 

INTEGER  I, I I, J, K,NPNTS, RES, METHOD, STARTND, STOPND, MODE, NRES, LB I AS 

INTEGER  GBIAS 

CHARACTER''- 1  CHAR ,  CHARl ,  CHAR2  ,  CHARS  ,  CHAR4 ,  CHARS 

CHARACTER-'n2  NAME 

COMMON/BLKl/MESH 
C0MM0N/BLK9/CHAR2,  CHARS,  CHAR4 


C  =  2. 99792SE+08 

PI  =  4. 

0''^ATAN(1.  0) 

READ(S 

1070)  NAME 

READ(3 

1000)  CHAR 

READ(S 

1000)  CHAR2 

READ(S 

1000)  CHARS 

READ(S 

1000)  CHAR4 

READ(3 

1000)  CHARS 

RE AD (3 

^0  MRESW 

READ(3 

'>)  DISTW 

RE AD (3 

^0  DPER 

RE AD (3 

''0  LB  IAS 

READ(S 

■^0    GBIAS 

READ(S 

''0  MAXD 

READ(S 

'■n    RES 

RE AD (3 

■>'0   NRES 

READ( 3 

,-^n    METHOD 

READ( 3 

,-'0  STARTND 

RE AD (3 

,^0  STOPND 

READ(3 

,■>'')    XORIGIN 

RE AD (3 

,^0  YORIGIN 

RE AD (3 

,''0  DD 

RE AD (3 

,")  A 
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READ(3,''0  B 

ALPHA  =  CMPLX(A,B) 

READ(3,'^)  A 

READO,''^)  B 

BETA  =  CMPLX(A,B) 

READ(3,1000)  CHARl 

IF(METHOD.EQ.  1)  THEN 

WRITE(6,''0  'METHOD  1  SELECTED  -  OBJECT  BISECTION  ' 
ELSEIF(METHOD. EQ. 2)  THEN 

WRITE(6,*)  'METHOD  2  SELECTED  -  CONNECTED  NODES  ' 
ELSEIF(METHOD.EQ. 3)  THEN 

WRITE (6,''-)  'METHOD  3  SELECTED  -  SELECTED  STOP  NODE  ' 
ELSE IF( METHOD.  EQ.  4)  THEN 

WRITE (6,*)  'METHOD  4  SELECTED  -  CONNECTED/ SELECTED  STOP  NODE' 
ELSEIF(METHOD.EQ. 5)  THEN 

WRITE(6,*)  'METHOD  5  SELECTED  -  EQUALLY  SPACED  CONNECTD  NODE' 
ELSE 

WRITE(6,*)  'METHOD  6  SELECTED  -  SELECTED  START  AND  STOP  NODE' 
END  IF 

IF((CHAR2.EQ.  'l').OR.  (CHAR2.EQ.  'i'))  THEN 

WRITE(6,*)  'INTERMEDIATE  MATRIX  FILE  GENERATION   -  ENABLED' 
ELSE 

WRITE(6,*)  'INTERMEDIATE  MATRIX  FILE  GENERATION   -  DISABLED' 

ENDIF 

IF(  ( CHARS.  EQ.  'D').0R.  (CHARS.  EQ.  'd'))  THEN 

WRITE  ( 6, '"0  'DISSPLA  FORTRAN  PROGRAM  GENERATION   -  ENABLED' 
ELSE 

WRITE(6,'V)  'DISSPLA  FORTRAN  PROGRAM  GENERATION   -  DISABLED' 

ENDIF 

IF((CHAR4.EQ.  'U').OR.  (CHAR4.EQ.  'u'))  THEN 

WRITE(6,''0  'UNIFORM  MATERIAL  SPECIFIED  (NO  INTERFACE)  ' 
ELSE 

WRITE(6,''0  'MATERIAL  SPECIFIED  WITH  A  VACUUM  AROUND  OBJECT' 
ENDIF 


IF((CHAR5.EQ.  •M').0R.  (CHARS.  EQ.  'm'))  THEN 

WRITE(6,''0  'MESH  GENERATION  <«  ONLY  >» 

ELSE 

WRITE ( 6, 'V)  'MESH  GENERATION  AND  FE  PROGRAM 

ENDIF 


ENABLED' 
ENABLED' 


IF(  (CHARl.  EQ.  'P').OR.  (CHARl.  EQ.  'p'))  THEN 

READ(3,''^)  EO 

READ(S,--)  FO 

LAMBDA  =  C/FO 

KO  =  2 "PI /LAMBDA 

WRITE(6,*)  'PLANEWAVE  BOUNDARY  VALUE  GENERATION   -  ENABLED' 

WRITE(6,^>)  'AMPLITUDE(EO)  =  ',  EO,',  WAVENUMBER(KO)  =  '  ,K0 

WRITE (6,*) 'WAVELENGTH    =  ', LAMBDA,',  FREQUENCY(FO)   =  ' ,F0 
ELSEIF(  (CHARl.  EQ.  'C').OR.  (CHARl.  EQ.  'c'))  THEN 

READ(3,--0  EO 

RE  AD  (3,^-)  MODE 
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READ(3,*)  FO 
LAMBDA  =  C/FO 
KO  =  2''^PI'-LAMBDA 

VRITE(6,''0  'CLYINRICAL  BOUNDARY  VALUE  GENERATION   -  ENABLED' 
WRITE(6,''0  'AMPLITUDE(EO)  =  ',  EO,',  MODE  NUMBER  =  ',  MODE 
WR I  TE(  6, ''0 'WAVELENGTH    =  ', LAMBDA,',  FREQUENCY(FO)   =  '  ,F0 
ELSE 

WRITE(6,''f)  'ALL  BOUNDARY  VALUE  GENERATION  METHODS  -  DISABLED' 
DO  5,  K  =  1,  RES 

READ(3,''0  BCOND(K) 
5  CONTINUE 

END  IF 
C 

C     POLAR  COORDINATE  INPUT  ROUTINE 
C 

IF((CHAR.EQ.  'P').OR.  (CHAR.  EQ.  'p'))  THEN 
READ(3,1020)  THETA,  RADIUS 
SFAC  =  2*PI^^DD/RADIUS 

MESH(0,4)  =  (RADIUS^'>-SIN(THETA''^PI/180.0))''^SFAC  +  XORIGIN 
MESH(0,5)  =  (RADIUS'>COS(THETA^'^PI/180.  0))''^SFAC  +  YORIGIN 
DO  10,  J  =  1,  RES-1 

READ(3,1020)  THETA,  RADIUS 

MESH(J,4)  =  (RADIUS->SIN(THETA''^PI/180.  0))*SFAC  +  XORIGIN 
MESH(J,5)  =  (RADIUS''^C0S(THETA''^PI/180.  0))''^SFAC  +  YORIGIN 
10         CONTINUE 

WRITE(6,'V)  'POLAR  COORDINATE  INPUT  SELECTED  ' 
C 

C     RECTANGULAR  COORDINATE  INPUT  ROUTINE 
C 

ELSEIF((CHAR.EQ.  'R').0R.  (CHAR.  EQ.  'r'))  THEN 
READ(3,1010)  MESH(0,4),  MESH(0,5) 
SFAC  =  2"PI"DD/((MESH(0,4)^--^--2+MESH(0,5)''-'''2)  — 0.  5) 
MESH(0,4)  =  MESH(0,4)-''^SFAC  +  XORIGIN 
MESH(0,5)  =  MESH(0,5)''^SFAC  +  YORIGIN 
DO  20,  J  =  1,  RES-1 

READ(3,1010)  MESH(J,4),  MESH(J,5) 
MESH(J,4)  =  MESH(J,4)'VSFAC  +  XORIGIN 
MESH(J,5)  =  MESH(J,5)^^SFAC  +  YORIGIN 
20         CONTINUE 

WRITE(6,^0  'RECTANGULAR  COORDINATE  INPUT  SELECTED  ' 
ELSE 

WRITE(6,''')  'INPUT  DATA  FILE  COORDINATE  SPECIFICATION  ERROR' 
ENDIF 

WRITE ('-,1100)  DD,  SFAC 
MESH(RES,4)  =  MESH(0,4) 
MESH(RES,5)  =  MESH(0,5) 
NPNTS  =  J 

WRITE(6,''0  'NODE  AT  0  DEGREES  (X/Y  COORDINATES)  =  ',  MESH(0,4), 
CMESH(0,5) 
WRITE(6,^-)  'X,  Y  OFFSETS  =  ',  XORIGIN,  YORIGIN 
WRITE(6,1090)  ALPHA,  BETA 
MRES  =  MRESW^>2^'^PI 

WRITE(6,'>)  'MESH  RESOLUTION       =  ',MRESW,'  WAVELENGTHS' 
DIST  =  DISTW^'-2^'--PI 

WRITE(6,''-)  'CONTOUR  DISTANCE      =  ',DISTW,'  WAVELENGTHS' 
WRITE(6,-0  'NUMBER  OF  DATA  POINTS  =  ' ,RES 
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IF(  ( METHOD.  EQ.  3).  OR.  (METHOD.  EQ.  6)  )  THEN 


END  IF 

WRITE(4 
WRITE(4 

C    55 
WRITE(4 

C    80 
WRITE (4 

C    90 
WRITE (4 

C  100 
WRITE (4 

C  110 
WRITE(4 

C  120 
WRITE (4 

C  130 
WRITE (4 

C  140 
WRITE (4 

C  150 
WRITE (4 

C  160 
WRITE (4 

C  170 
WRITE (4 

C  180 
WRITE (4 

C  190 
WRITE (4 

C  200 
WRITE (4 

C  210 
WRITE(4 

C  220 
WRITE (4 

C  230 
WRITE (4 
WRITE (4 
WRITE(4 
WRITE(4 
WRITE (4 
WRITE(4 
WRITE (4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE (4 
WRITE (4 
WRITE(4 
WRITE(4 
WRITE (4 


WRITE (6,'"^)  'START  NODE 
WRITE (6,'^)  'STOP  NODE 


,v)  '17' 
■>)  '10 

,V)  '2 
,V)  '2 

*)  '2 

,v)  '2 

,v)  '2 

,v)  '2 

,v)  -2 

-,v)  '2 

,v)  '2 

,v)  '2 

*)    '2 
,v)  '2 

*)    '2 

,v)  '2 

,v)  '2 

,v)  '2 

2' 
1070)  NAME 

''")  'Mesh  Resolution' 
V.)  -L' 
1030)  MRESW 

''0  'Contour  Distance' 

1030)  DISTW 

^)  'L' 

")  'Number  of  Points' 

^v)  'L' 

1040)  RES 

V.)  'L' 

")  'Method' 
^0  'L' 


1' 
1' 
2' 
1* 
2' 
1' 
2' 
1' 
2' 
1' 
2' 
1' 
2' 
1' 
2' 
1' 


=  ' , STARTND 

=  ' , STOPND 

0 

130 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 

0 

365 
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30 


C 

1000 
1010 
1020 
1030 
1040 
1050 
1060 
1070 
1080 
1090 

1100 


WRITE (4 
WRITE (4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE (4 
WRITE (4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE(4 
WRITE (4 
WRITE (4 
WRITE (4 
WRITE (4 
WRITE(4 


,1040)  Method 

,''0  'X  Origin' 

,'0  'L' 

,1050)  XORIGIN 

,*)    'L' 

,-'0  'Y  Origin' 

,1050)  YORIGIN 
,'>)  'L' 
,■''0    'Alpha' 
,''-)  'L' 
,1060)  ALPHA 
,'-)  'L' 
,''0  'Beta' 
,^v)  'L' 
,1060)  BETA 
,'■)  'L' 
,''0  '9' 


DO  30,  II  =  0,  NPNTS-1 

WRITE(2,'^)  MESH(II,4),  MESHCII,5) 
CONTINUE 

WRITE (2,---)  MESH(0,4),  MESH(0,5) 
WRITE(2,*)  '999992,  999991  ' 
WRITE(2,-''-)  '999990,  999990  ' 

FORMAT(Al) 

FORMAT(10X,F12.  8,4X,F12.  8) 
F0RMAT(3X,I3,5X,F12.  8) 
F0RMAT(F8.  6) 
F0RMAT(I4) 
F0RMAT(1X,F8.  6) 
F0RMAT(1X,F8.  6,'  -  J' ,F8.  6) 
F0RMAT(A12) 
FORNATC/) 

FORMAT( IX, 'ALPHA  =  ',F8.4,2X,'+  J  ',F6. 4,' 
C  '  ,F6.4) 
FORMATC IX, 'DESIRED  DISTANCE  =  ',F8.5,', 
RETURN 
END 


BETA  =  'F8.4,2X,'+  J 
SCALE  FACTOR  =  ' ,F8.5) 


SUBROUTINE  ROTATE ( NPNTS , METHOD , STARTND , STOPND) 

THIS  SUBROUTINE  ROTATES  THE  SURFACE  POINTS  TO  ALLOW  FOR 
A  REARRANGING  OF  THE  START  NODE  (FIRST  DATA  POINT). 

REAL  MESH(0: 200,5) 

INTEGER  I,  NPNTS,  METHOD,  STARTND,  STOPND 

COMMON/BLKl/MESH 


IFCMETHOD.EQ.  6)  THEN 

DO  10,  I  =  0,  NPNTS-1 

MESH(I,1)  =  MESH(I,4) 
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MESH(I,2)  =  MESH(I,5) 
10         CONTINUE 
C 

DO  20,  I  =  STARTND-1,  NPNTS-1 

MESH(I-(STARTND-1),4)  =  MESH(I,1) 
MESH( I -(STARTND-1), 5)  =  MESH(I,2) 
20         CONTINUE 
C 

DO  30,  1=1,  STARTND 

MESH(I+NPNTS-STARTND,4)  =  MESH( 1-1,1) 
MESH(I+NPNTS-STARTND,5)  =  MESH( 1-1,2) 
30         CONTINUE 
ENDIF 

STOPND  =  STOPND  -  STARTND  +  1 
IF(STOPND. GT.  NPNTS)  THEN 

STOPND  =  STOPND  -  NPNTS 
ELSEIF(STOPND. LT. 0)  THEN 

STOPND  =  NPNTS  +  STOPND  +1 
ELSEIF( STOPND. EQ. 0)  THEN 

STOPND  =  1 
ELSE 

STOPND  =  STOPND 
ENDIF 
RETURN 
END 
C 
C 

SUBROUTINE  BOUND( MRES , NPNTS , METHOD , STOPND ) 
C 

C     THIS  SUBROUTINE  REORDERS  THE  THE  SURFACE  POINTS 
C     BASED  ON  THE  DESIRED  INPUT  MESH  RESOLU  INN.  THE 
C     OBJECT  IS  ALSO  BISECTED. 
C 

REAL  PERIM,  DIST,  MESH( 0: 200 ,5) ,  TEMP,  TEMPI,  MRES,  MRESN,  MRESLW 
REAL  MRESNL,  MRESNR,  TEMPR,  TEMPL,  DISTB,  DZ,  PI,  MRESW,  MRESRW 
INTEGER   I,  PERND,  BIND,  NPNTS , STARTND ,  STOPND,  METHOD,  NEWND,  J 
INTEGER  M,  L 
COMMON/BLKl/MESH 
C0MM0N/BLK2/PERND,BIND 
C 

PI  =  4.  0-''-ATAN(l.  0) 
PERIM  =  0.  0 
C 

DO  10,  I  =  0,  NPNTS-1 

DIST  =  SQRT((MESH(I+1,4)-MESH(I,4))**2+(MESH(I+1,5)  - 
C     MESH(I,5))^--''-2) 

PERIM  =  PERIM  +  DIST 
MESH( 1+1,3)  =  PERIM 
10    CONTINUE 

WRITE(6,''0  'PERIMETER  LENGTH      =  ',  PERIM 

PERND  =  NINT(PERIM/MRES  -  AMOD(PERIM/MRES , 2.  0) ) 

BIND  =  (PERND  -  2)/2 

WRITE(6,''0  'PERIMETER  NODE  //      =  ',  PERND 

MRESW  =  MRES/(2-''-PI) 

WRITE(6,-'0  'YOUR  REQUESTED  MESH  RESOLUTION  OF  ', MRESW 

IF(  ( METHOD.  EQ.  3).  OR.  (METHOD.  EQ.  4).  OR.  (METHOD.  EQ.  5).  OR. 
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C(METHOD.EQ.  6))  THEN 
C 

MRESNL  =  (PERIM  -MESH(ST0PND-1 ,3) )/FL0AT(PERND/2) 
MRESLW  =  MRESNL/ (2--'-P  I) 

MRESNR  =  MESH(ST0PND-l,3)/FL0AT(PERND/2) 
MRESRW  =  MRESNR/ (2-PI) 

WRITE(6,*)  '.  .  .   HAS  BEEN  MODIFIED  TO  ...  ' 
WRITE(6,''>-)  'LEFT  SIDE  OF  SEGMENT  .  .  .  ', MRESLW 
WRITE(6,'^)  'RIGHT  SIDE  OF  SEGMENT  .  .  ,  ', MRESRW 
ELSE 

MRESN  =  PERIM/FLOAT(PERND) 
MRESW  =  MRESN/(2'-PI) 

WRITE(6,*)  '.  .  .   HAS  BEEN  MODIFIED  TO  .  .  .  ', MRESW 
ENDIF 
C 

C     PERIMETER  NODE  INITIALIZATION  (X,Y  COORD) 
C 

IF( ( METHOD.  EQ.  1).0R.  (METHOD.  EQ.  2)  )  THEN 
DO  30,  1=0,  PERND-1 
TEMP  =  MRESN" I 
J  =  1 
20  IF(TEMP.GT. MESH(J,3))  THEN 

J  =  J  +  1 
GOTO  20 
ENDIF 
TEMPI  =  (TEMP  -  MESH(J-1,3))/(SQRT((MESH(J,4)-MESH(J-1,4)) 
C     *''-2  +  (MESH(J,5)  -  MESH(J-1,5))''"'^2)) 

MESH( 1+1,1)  =  MESH(J-1,4)  +  TEMP1"(MESH( J,4) -MESH( J-1 ,4) ) 
MESH(I+1,2)  =  MESH(J-1,5)  +  TEMP1'KMESH(  J,5) -MESH(  J-1 ,5) ) 
30         CONTINUE 
C 

ELSE 

DO  60,  1=0,  PERND-1 
TEMPR  =  MRESNR^-- 1 

TEMPL  =  PERIM  -  MRESNL''-(PERND  -  I) 
IF( TEMPR. LE. MESH( STOPND- 1 , 3 ) )  THEN 
J  =  1 
40  IF(TEMPR. GT.  MESH(J,3))  THEN 

J  =  J  +  1 
GOTO  40 
ENDIF 
TEMPI  =  (TEMPR  -  MESH( J-1 , 3) )/( SQRT( (MESH( J,4) -MESH( J-1 ,4) ) 
C     *,V2  +  (MESH(J,5)  -  MESH(J-1,5))''"'^2)) 

MESH(I+1,1)  =  MESH(J-1,4)  +  TEMP1^'-(MESH(  J,4) -MESH(  J-1 ,4)  ) 
MESH(I+1,2)  =  MESH(J-1,5)  +  TEMP1''-(MESH(  J,5) -MESH(  J-1 ,5) ) 
ELSE 

J  =  1 
50  IF(TEMPL.  GT.  MESH(J,3))  THEN 

J  =  J  +  1 
GOTO  50 
ENDIF 
TEMPI  =  (TEMPL  -  MESH( J-1 ,3) )/(SQRT( (MESH( J,4)-MESH( J-1,4) ) 
C     '•-*2  +  (MESH(J,5)  -  MESH(J-l,5))''-"2)) 

MESH(I+1,1)  =  MESH(J-1,4)  +  TEMP1'-(MESH(  J,4)-MESH(  J-1 ,4) ) 
MESH( 1+1,2)  =  MESH(J-1,5)  +  TEMP1"(MESH( J,5)-MESH( J-1 ,5) ) 
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ENDIF 
60         CONTINUE 

ENDIF 
C 
C 

C     BISECTION  NODE  INITIALIZATION  (X,Y  COORD) 
C 

WRITE(6,*)  'BISECTION  NODE  #      =  ' ,BIND 
IF(  (METHOD.  EQ.  1).0R.  (METHOD.  EQ.  3)  )  THEN 
DO  70,  1=1,  BIND 

TEMPI  =  I/FLOAT(BIND  +  1) 

MESH(I+PERND,1)  =  MESH(1,1)  +  TEMPl*(MESH(BIND+2, 1)  - 
CMESH(1,1)) 

MESH(I+PERND,2)  =  MESH(1,2)  +  TEMPl'V(MESH(BIND+2,2)  - 
CMESH(1,2)) 
70         CONTINUE 
C 

ELSE 

DO  80,  1=2,  PERND+1 

NESH(PERND+I-1,1)  =  MESH(PERND+2-1 , 1)  +  0.  5''^(MESH(  1 , 1) 
C  MESH(PERND+2-I,l)) 

MESH(PERND+I-1,2)  =  MESH(PERND+2-1 ,2)  +  0.  5''-(MESH(  1 ,2) 
C  MESH(PERND+2-I,2)) 

80         CONTINUE 
C 
C 

ENDIF 
C 

IF(  ( METHOD.  EQ.  5 ) .  OR.  ( METHOD.  EQ.  6 ) )  THEN 
DO  90,  M  =  1,  BIND 

MESH(M,4)  =  MESH(PERND+M,1) 
MESH(M,5)  =  MESH(PERND+M,2) 
90         CONTINUE 

DISTB  =  0.0 
MESH(0,3)  =  0.  0 
DO  100,  L  =  1,  BIND+1 
IF(L.  EQ.  1)  THEN 

DIST  =  SQRT((MESH(1,1)  -  MESH(  1 ,4)  )^'^^>-2  + 
C(MESH(1,2)  -  MESH(l,5))->--"^2) 

ELSEIF(L.EQ.  BIND+1)  THEN 

DIST  =  SQRT((MESH(BIND,4)  -  MESH(BIND+2 , 1) ) 
C-VV.-2  +  (MESH(BIND,5)  -  MESH(BIND+2 ,2)  )^'-''^2) 
ELSE 

DIST  =  SQRT((MESH(L-1,4)  -  MESH(L,4)) 
C^v,V2  +  (MESH(L-1,5)  -  MESH(L,5)  )**2) 
ENDIF 

DISTB  =  DISTB  +  DIST 
MESH(L,3)  =  DISTB 
100        CONTINUE 

DZ  =  DISTB/(BIND  +  1.  0) 
WRITE(6,->)  'DZ  SPACING  =  '  ,  DZ 

C 

DO  120,  I  =  1,  BIND 
TEMP  =  DZ^^I 
J  =  1 
110  IF(TEMP.  GT. MESH(J,3))  THEN 
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J  =  J  +  1 

GOTO  110 
ENDIF 

IF(J.EQ.  1)  THEN 
TEMPI  =  TEMP/(SQRT((MESH(1,1)  -  MESH( 1 ,4) )**2  + 
C     (MESH(1,2)  -  MESH(1,5))^"'^2)) 

MESH(PERND+I,1)  =  MESH(1,1)  +  TEMP 1*(MESH( 1,4) 
C     -  MESH(1,1)) 

MESH(PERND+I,2)  =  MESH(1,2)  +  TEMP1*(MESH(1 ,5) 
C     -  MESH(1,2)) 
ELSE 
TEMPI  =  (TEMP  -  MESH(J-1,3))/(SQRT((MESH(J,4)  - 
C     MESH(J-l,4))''"'-2  +  (MESH(J,5)  -  MESH(  J-1 ,5)  )''"'f2) ) 
MESH(PERND+I,1)  =  MESH(J-1,4)  +  TEMP1*(MESH( J,4) 
C     -  MESH(J-1,4)) 

MESH(PERND+I,2)  =  MESH(J-1,5)  +  TEMP1*(MESH( J,5) 
C     -  MESH(J-1,5)) 
ENDIF 
120        CONTINUE 
ENDIF 
RETURN 
END 
C 
C 

SUBROUTINE  NORMAL 
C 

C     THIS  ROUTINE  COMPUTES  THE  X  AND  Y  COMPONENTS 
C     OF  THE  OUTWARD  UNIT  NORMAL  AT  EACH  SURFACE  POINT. 
C 

REAL  DR,  DZ,  DL,  MESH(0:  200 ,5) 
INTEGER  I,  PERND,  BIND 
COMMON/BLKl/MESH 
C0MM0N/BLK2/PERND , BIND 
DO  10,  I  =  1,  PERND 
IFCI.EQ.  1)  THEN 

DR  =  MESH(2,1)    -   MESH( PERND, 1) 
DZ   =  MESH(2,2)    -   MESH(PERND, 2) 
DL  =   SQRT(DR--'^DR+DZ"DZ) 
MESH(1,3)=-DZ/DL 
MESH(1,4)=DR/DL 
ELSEIFd.EQ.  PERND)   THEN 

DR  =  MESH(1,1)  -  MESH( PERND- 1,1) 
DZ  =  MESH(1,2)  -  MESH( PERND- 1,2) 
DL  =  SQRT(DR--'-DR+DZ^'-DZ) 
MESH( PERND, 3)=-DZ/DL 
MESH( PERND, 4)=DR/DL 


ELSE 


ENDIF 
10    CONTINUE 
RETURN 


DR  =  MESH( 1+1,1)  -  MESH( 1-1,1) 
DZ  =  MESH(I+1,2)-  MESH(I-1,2) 
DL  =  SQRT(DR"DR+DZ^^DZ) 
MESH(I,3)=-DZ/DL 
MESH(I,4)=DR/DL 
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END 
C 
C 

SUBROUTINE  NODSET( METHOD , NABC , NBMX , UNK) 
C 

C     THIS  ROUTINE  USES  A  TWO-SWEEP  TECHNIQUE  TO  COMPUTE  THE 
C     NUMBER  OF  NODES  ALONG  EACH  NODE  ROW,  I.  A  MESH 
C     ORIENTATION  ATTRIBUTE  IS  ALSO  SET. 
C 

C     SET  Z-AXIS  SPACING  AND  ENDPOINT  NODES 
C 

REAL  DZ,  ZZ,  MESH(0: 200,5),  D,  DIST,  DISTB 

INTEGER  I,  J,  NODES(200),  M0R(200),  PERND,  NOLD,  NNEW,  BIND 
INTEGER  L,  METHOD,  NABC( 100,3),  NBMX,  UNK 
CHARACTER'n  CHAR2,  CHARS,  CHAR4 
COMMON/BLKl/MESH 
C0MM0N/BLK2/PERND,BIND 
C0MM0N/BLK3/N0DES 
C0MM0N/BLK6/M0R 
COMMON/ BLK9/CHAR2, CHARS, CHAR4 
C 

UNK  =  0 

NBMX  =  -99 

IF(  (METHOD.  EQ.  1).0R.  (METHOD.  EQ.  3))  THEN 

DZ  =  (SQRT((MESH(1,2)  -  MESH(BIND+2 ,2) )**£  + 
C     (MESH(1,1)  -  MESH(BIND+2,1))''"'^2))/(BIND  +  1.0) 
ELSE 

DISTB  =  0.0 
DO  10,  L  =  1,  BIND+1 
IF(L.  EQ.  1)  THEN 

DIST  =  SQRT((MESH(1,1)  -  MESH(PERND+1 , 1) )**2  + 
C(MESH(1,2)  -  MESH(PERND+l,2))"^>-2) 
ELSEIF(L.EQ.  BIND+1)  THEN 

DIST  =  SQRT((MESH(PERND+BIND,1)  -  MESH(BIND+2 , 1) ) 
CVr,V2  +  (MESH(PERND+BIND,2)  -  MESH(BIND+2 ,2)  )'>''^2) 
ELSE 

DIST  =  SQRT((MESH(PERND+I-1,1)  -  MESH(PERND+I , 1) ) 
Cv.-vr2  +  (MESH(PERND+I-1,2)  -  MESH(  PERND+1 ,2)  )^''''2) 
ENDIF 

DISTB  =  DISTB  +  DIST 
10         CONTINUE 

DZ  =  DISTB/(BIND  +  1.0) 
ENDIF 

WRITE(6,''0  'BISECTION  SPACING     =  '  ,  DZ 
NODES(l)  =  2 
N0DES(2)  =  S 
N0DES(BIND+2)  =  2 
N0DES(PERND+1)  =  2 
NODES( PERND)  =  3 
C     PERFORMING  FORWARD-SWEEP 
DO  20  I  =  3,  BIND+1 

NOLD  =  NODES(I-l) 

D  =  SQRT((MESH(I,1)  -  MESH( I+PERND-1 , 1) )**2  + 
C     (MESH(I,2)-MESH(I+PERND-l,2))^^''-2) 
NNEW  =  INT(0. 1  +  D/DZ)  +  2 
NODES(I)  =  NNEW 


IF  (NNEW. GT. NOLD+1)  NODES(I)  =  NOLD  +  1 
IF  (NNEW.  LT. NOLD- 1)  NODES(I)  =  NOLD  -  1 
IF  (NODES(I).LT.  3)  NODES(I)  =  3 
20    CONTINUE 
C 

J  =  PERND  +  2 

DO  30,  I  =  PERND- 1,  BIND+3,  -1 
NOLD  =  N0DES(I+1) 

D  =  SQRT((MESH(I,1)  -  MESH(  J,  1)  )*''f2  + 
C     (MESH(I,2)-MESH(J,2))''"'-2) 
NNEW  =  INT(0. 1  +  D/DZ)  +  2 
NODES(I)  =  NNEW 

IF  (NNEW.  GT. NOLD+1)  NODES(I)  =  NOLD  +  1 
IF  (NNEW. LT. NOLD- 1)  NODES(I)  =  NOLD  -  1 
IF  (NODES(I).LT.  3)  NODES(I)  =  3 
J  =  J  +  1 
30    CONTINUE 
C 

C     BACK-SWEEP  TO  RESET  LAST  NODES  IF  NEEDED 
C 

I  =  BIND  +  2 
40    1=1-1 

IF  (I.EQ.  2)  GO  TO  50 

IF  (N0DES(I).LE.N0DES(I+1)+1)  GO  TO  60 
NODES(I)  =  N0DES(I+1)  +  1 
GO  TO  40 
50    CONTINUE 

WRITE(6,^0  '  PROGRAM  ABORTED  IN  NODSET  RIGHT  SIDE  BACKSWEEP  ' 
STOP 
60    CONTINUE 
C 

I  =  BIND  +  2 
70    1=1+1 

IF  (I.EQ.  PERND)  GO  TO  80 
IF  (N0DES(I).LE.N0DES(I-1)  +  1)  GO  TO  90 
NODES(I)  =  NODES(I-l)  +  1 
GO  TO  70 
80    CONTINUE 

WRITE(6,*)  '  PROGRAM  ABORTED  IN  NODSET  LEFT  HALF  BACKSWEEP  ' 
STOP 
90    CONTINUE 
C 

C     FORWARD  SWEEP  TO  LOAD  MESH  ORIENTATION  ARRAY,  MOR 
C 

MOR(l)  =  0 
M0R(BIND+2)  =  0 
M0R(PERND+1)  =  0 
C 

DO  100,  1=2,  BIND+1 

IF(N0DES(I+1).GT. NODES(I))  THEN 

MOR(I)  =  0 
ELSEIF(N0DES(I+1).LT.  NODES(I))  THEN 

MOR(I)  =  1 
ELSE 

IF(N0DES(I+2).GT. NODES(I))  THEN 
MOR(I)  =  0 
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ELSEIF(N0DES(I+2).LT. NODES(I))  THEN 

MOR(I)  =  1 
ELSE 

MOR(I)  =  MOR(I-l) 
ENDIF 
ENDIF 
100   CONTINUE 
C 

DO  110,  I  =  PERND,  BIND+3,  -1 

IF(N0DES(I-1).GT. NODES(I))  THEN 

MOR(I)  =  0 
ELSEIF(N0DES(I-1).LT.  NODES(I))  THEN 

MOR(I)  =  1 
ELSE 

IF(N0DES(I-2).GT.  NODES(I))  THEN 

MOR(I)  =  0 
ELSEIFC NODES( I -2 ) . LT. NODES( I ) )  THEN 

MOR(I)  =  1 
ELSE 

MOR(I)  =  M0R(I+1) 
ENDIF 
ENDIF 
110   CONTINUE 
C     LOAD  NABC  ARRAY 

DO  120,  I  =  1,  BIND+2 
IF(I.EQ.  1)  THEN 

NABC(1,1)  =  0 

NABC(1,2)  =  1 

NABC(1,3)  =  3 

ELSEIF(I.LE.BIND)  THEN 

NABC(I,1)  =  NABC(I-1,2) 
NABC(I,2)  =  NABC(I-1,3) 

NABC(I,3)  =  N0DES(PERND-I+1)  +  N0DES(I+1)  -  3 
ELSEIFC  I.  EQ.BIND+1)  THEN 

NABC(I,1)  =  NABC(I-1,2) 
NABC(I,2)  =  NABC(I-1,3) 
NABC(I,3)  =  1 
ELSE 

NABC(I,1)  =  NABC(I-1,2) 
NABC(I,2)  =  NABC(I-1,3) 
NABC(I,3)  =  0 
ENDIF 
120   CONTINUE 
C 

DO  130,  I  =  1,  BIND+2 

IF((CHAR2.EQ.  •l').OR.  (CHAR2.EQ.  'i'))  THEN 

WRITE(  12 ,''-)  I  ,NABC(  1 , 1)  ,NABC(  1,2)  ,NABC(  1,3) 
ENDIF 

UNK  =  UNK  +  NABC(I,2) 
IF( NABC ( 1,2). GE.NBMX)  THEN 

NBMX  =  NABC (1,2) 
ENDIF 
130   CONTINUE 
C 

WRITE('^-'^)  'MAXIMUM  UNKNOWN  WIDTH  =  '  ,NBMX 
WRITEC^'S^'O  'TOTAL  #  OF  UNKNOWNS   =  '  ,UNK 
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WRITE(''s*)  '  ' 

RETURN 

END 
C 
C 

c 

SUBROUTINE  SORTER( I ,LEL,LAND) 
C 

C     THIS  SUBROUTINE  GENERATES  A  MESH  ROW  FOR  THE  INPUT  ROW  I. 
C     LOADING  OF  THE  LOCAL  NODE-ELEMENT  CONNECTION  MATRICES  LND  AND 
C     NDL  FOR  ELEMENTS  BETWEEN  I  AND  I+l  NODE  ROWS.  REFERENCE  IS  TO 
C     THE  LEFT  SIDE  OF  THE  Ith  ROW  OR  VECTOR. 
C 

INTEGER  NODES( 200),  MOR(200),  LND(0: 200,3) ,  NDL(200,4),  NCT(200) 

INTEGER  I,  PERND,  BIND,  LEL,  LAND,  J,  K,  LL,  JJ,  NN,  KK,  N,  NX,  N2 

INTEGER  NDMX,  NDSIL,  NDS2L,  NDSIR,  NDS2R,  LMX 

C0MM0N/BLK2/PERND , BIND 

C0MM0N/BLK3/N0DES 

C0MM0N/BLK4/LND,NDL,NCT 

COMMON/ BLK6/M0R 

IF(I.EQ.  BIND+2)  THEN 

WRITE (6,''^)  'ERRORED  OUT  IN  SUBROUTINE  SORTER,  YOU  ATTEMPTED' 
WRITE(6,*)  '      TO  CALL  SORTER  WITH  I  =  BIND  +  2!  ' 
WRITE(6,'^)  '  THIS  ROW  HAS  NO  ELEMENTS' 

RETURN 
END  IF 
DO  20,  J  =  0,  200 

DO  10,  K  =  1,  3 

LND(J,K)  =  0 
10         CONTINUE 
20    CONTINUE 

NDMX  =  200 

NDSIL  =  N0DES(PERND+2-I) 
NDS2L  =  N0DES(PERND+1-I) 
NDSIR  =  NODES(I) 
NDS2R  =  N0DES(I+1) 

LMX  =  NDSIL  +  NDSIR  +  NDS2L  +  NDS2R  -  2 
C 
C 

C     LEFTSIDE  OF  THE  BISECTION  SEGMENT 
C 

C     TOP  ROW 
C 

IFCI.EQ.  1)  THEN 

LND(1,1)  =  4 
LND(1,2)  =  1 
LND(1,3)  =  3 
LND(2,1)  =  1 
LND(2,2)  =  4 
LND(2,3)  =  2 
LND(3,1)  =  5 
LND(3,2)  =  2 
LND(3,3)  =  4 
LND(4,1)  =  5 
LND(4,2)  =  2 
LND(4,3)  =  6 
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LND(5,1)  =  1 

LND(5,2)  =  6 

LND(5,3)  =  2 

LND(6,1)  =  6 

LND(6,2)  =  1 

LND(6,3)  =  7 

LEL  =  6 

LAND  =  7 

c 

c 

BOTTOM  ROW 

c 

ELSEIFd.EQ.  BIND+1 

LND(1,1)  =  2 

LND(1,2)  =  6 

LND(1,3)  =  1 

LND(2,1)  =  6 

LND(2,2)  =  2 

LND(2,3)  =  7 

LND(3,1)  =  3 

LND(3,2)  =  7 

LND(3,3)  =  2 

LND(4,1)  =  3 

LND(4,2)  =  7 

LND(4,3)  =  4 

LND(5,1)  =  6 

LND(5,2)  =  4 

LND(5,3)  =  7 

LND(6,1)  =  4 

LND(6,2)  =  6 

LND(6,3)  =  5 

LEL  =  6 

LAND  =  7 

THEN 


EQUAL  NODE  NUMBERS 

ELSEIF(NDS1L. EQ.NDS2L)  THEN 
FOR  MOR  =  0  (LH  ORIENTATION) 

IF(M0R(PERND+2-I).EQ.  0)  THEN 


LND(1,1)  =  1 

LND(1,2)  =  NDSIL  +  NDSIR 

LND(1,3)  =  2 

LND(2,1)  =  NDSIL  +  NDSIR  +  1 

LND(2,2)  =  2 

LND(2,3)  =  NDSIL  +  NDSIR 

DO  40,  N  =  1,  NDSlL-2 

Nl  =  2-''N  +  1 

N2  =  Nl  +  1 

DO  30,  K  =  1,  3 

LND(N1,K)  =  LND(1,K)  +  N 

LND(N2,K)  =  LND(2,K)  +  N 

30 

CONTINUE 

40 

CONTINUE 

C 

FOR  MOR  = 
ELSE 

1  (RH  ORIENTATION) 

LND(1,1)  =  NDSIL  +  NDSIR 
LND(1,2)  =  1 
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LND(1,3)  =  NDSIL  +  NDSIR  +  1 

LND(2,1)  =  2 

LND(2,2)  =  NDSIL  +  NDSIR  +  1 

LND(2,3)  =  1 

DO  60,  N  =  1,  NDSlL-2 

Nl  =  2-->N  +  1 

N2  =  Nl  +  1 

DO  50,  K  =  1,  3 

LND(N1,K)  =  LND(1,K)  +  N 

LND(N2,K)  =  LND(2,K)  +  N 

50 

CONTINUE 

60 

CONTINUE 

ENDIF 

C 

c 
c 

LEFTHAND  MESH  ORIENTATION 

ELSEIF(M0R(PERND+2-I).EQ. 0)  THEN 

LND(1,1)  =  NDSIL  +  NDSIR  +  1 

LND(1,2)  =  1 

LND(1,3)  =  NDSIL  +  NDSIR 

LND(0,1)  =  0 

LND(0,2)  =  NDSIL  +  NDSIR 

LND(0,3)  =  1 

DO  80,  N  =  1,  NDSIL- 1 

Nl  =  2'>N 

N2  =  Nl  +  1 

DO  70,  K  =  1,  3 

LND(N1,K)  =  LND(0,K)  +  N 

LND(N2,K)  =  LND(1,K)  +  N 

70 

CONTINUE 

80 
C 
C 
C 

CONTINUE 

RIGHTHAND  MESH  ORIENTATION 

ELSE 

LND(1,1)  =  2 

LND(1,2)  =  NDSIL  +  NDSIR 

LND(1,3)  =  1 

LND(0,1)  =  NDSIL  +  NDSIR  -  1 

LND(0,2)  =  1 

LND(0,3)  =  NDSIL  +  NDSIR 

DO  100,  N  =  1,  NDSlL-2 

Nl  =  2''^N 

N2  =  Nl  +  1 

DO  90,  K  =  1,  3 

LND(N1,K)  =  LND(0,K)  +  N 

LND(N2,K)  =  LND(1,K)  +  N 

90 

CONTINUE 

100 

CONTINUE 

ENDIF 

IF((I.EQ.  1).0R.  (I.EQ.  BIND+1))  THEN 

GOTO  230 
ENDIF 
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LEL  =  N2 


RIGHTSIDE  OF  THE  BISECTION  SEGMENT 


EQUAL  NODE  NUMBERS 

IF(NDS1R.  EQ.  NDS2R)  THEN 
MOR  =  0  (LH  ORIENTATION) 
IF(MOR(I).EQ.  0)  THEN 


LND(LEL+1,1)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 

LND(LEL+1,2)  =  NDSIL 

LND(LEL+1,3)  =  NDSIL  +  NDSIR  +  NDS2L 

LND(LEL+2,1)  =  NDSIL  +  1 

LND(LEL+2,2)  =  NDSIL  +  NT)S1R  +  NDS2L 

LND(LEL+2,3)  =  NDSIL 

DO  120,  N  =  1,  NDSlR-2 

Nl  =  2'VN  +  1  +  LEL 

N2  =  Nl  +  1 

DO  110,  K  =  1, 

3 

LND(N1,K) 

=  LND(LEL+1,K)  + 

N 

LND(N2,K) 

=  LND(LEL+2,K)  + 

N 

110 

CONTINUE 

120 

c 

CONTINUE 
LEL  =  N2 
LAND  =  LMX 
MOR  =  1  (RH  ORIENTATION) 
ELSE 

LND(LEL+1,1)  =  NDSIL 

LND(LEL+1,2)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 

LND(LEL+1,3)  =  NDSIL  +  1 

LND(LEL+2,1)  =  NDSIL  +  NDSIR  +  NDS2L 

LND(LEL+2,2)  =  NDSIL  +  1 

LND(LEL+2,3)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 

DO  140,  N  =  1,  NDSlR-2 

Nl  =  2-'^N  +  1  +  LEL 

N2  =  Nl  +  1 

DO  130,  K  =  1, 

3 

LND(N1,K) 

=  LND(LEL+1,K)  + 

N 

LND(N2,K) 

=  LND(LEL+2,K)  + 

N 

130 

CONTINUE 

140 

CONTINUE     -, 
ENDIF 
LEL  =  N2 

C 

c 
c 

■  LAND  =  LMX 

LEFT  HAND  MESH  ORIENTATION 

ELSEIF(MOR(I).EQ.  0)  THEN 
LND(LEL+1,1)  =  NDSIL 
LND(LEL+1,2)  =  NDSIL 
LND(LEL+1,3)  =  NDSIL 
LND(LEL+2,1)  =  NDSIL 
LND(LEL+2,2)  =  NDSIL 


+  NDSIR  +  NDS2L  -  1 


NDSIR 

1 

NDSIR 


+  NDS2L 


+  NDS2L 
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LND(LEL+2,3)  =  NDSIL 
C 

DO  160,  N  =  1,  NDSlR-1 
Nl  =  2''-N  +  LEL  +  1 
DO  150,  K  =  1,  3 

LND(N1,K)  =  LND(LEL+1,K)  +  N 
150  CONTINUE 

160        CONTINUE 
C 

DO  180,  N  =  1,  NDSlR-2 
N2  =  2*N  +  LEL  +  2 
DO  170,  K  =  1,  3 

LND(N2,K)  =  LND(LEL+2,K)  +  N 


170 

CONTINUE 

180 

CONTINUE 

C 

LEL  =  Nl 

LAND  =  LMX 

C 

C 

RIGHT  HAND  MESH  ORIENTATION 

C 

ELSE 

LND(LEL+1,1)  =  NDSIL 

LND(LEL+1,2)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 
LND(LEL+1,3)  =  NDSIL  +  1 
LND(LEL+2,1)  =  NDSIL  +  NDSIR  +  NDS2L 
LND(LEL+2,2)  =  NDSIL  +  1 

LND(LEL+2,3)  =  NDSIL  +  NDSIR  +  NDS2L  -  1 
C 

DO  200,  N  =  1,  NDSlR-2 
Nl  =  2-'-N  +  LEL  +  1 
DO  190,  K  =  1,  3 

LND(N1,K)  =  LND(LEL+1,K)  +  N 
190  CONTINUE 

200        CONTINUE 
C 

DO  220,  N  =  1,  NDSlR-3 
N2  =  2'''N  +  LEL  +  2 
DO  210,  K  =  1,  3 

LND(N2,K)  =  LND(LEL+2,K)  +  N 
210  CONTINUE 

220        CONTINUE 
C 

LEL  =  Nl 
LAND  =  LMX 
C 
C 

ENDIF 
C 
C 
C 

C     ZEROING  NDL  AND  NCT  PRIOR  TO  FILL 
C 

230   DO  250,  N  =  1,  NDMX 
NCT(N)  =  0 
DO  240,  K  =  1,  4 
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240 

250 

C 

C 

C 


260 
270 
C 

C 
C 

C 
C 
C 
C 


NDL(N,K)  =  0 

CONTINUE 
CONTINUE 

SCANNING  LND  ARRAY  TO  LOAD  NDL 

DO  270  LL  =  1,  LEL 

DO  260  JJ  =  1,  3 

NN  =  LND(LL,JJ) 
NCT(NN)  =  NCT(NN)  +  1 
KK  =  NCT(NN) 
NDL(NN,KK)  =  LL 
CONTINUE 
CONTINUE 

END 


SUBROUTINE  FINDER( I ,LAND,DIST) 

THIS  SUBROUTINE  DETERMINES  THE  (X,Y)  COORDINATES  OF  EACH 
NODE  IN  THE  CALLING  Ith  ROW  OR  VECTOR  MESH. 


INTEGER  I,  J,  LAND,  PERND,  BIND 
REAL  MESH(0: 200,5),  NDP(200,2), 
COMMON/ ELK 1 /MESH 
C0MM0N/BLK2/PERND , BIND 
C0MM0N/BLK3/N0DES 
C0MM0N/BLK5/NDP 
IF(I.EQ.  1)  THEN 


,  NODES(200) 
DIST 


NDP(1 
NDP(1 
NDP(2 
NDP(2 
NDP(3 
NDP(3 
NDP(4 
NDP(4 
NDP(5 
NDP(5 
NDP(6 
NDP(6 
NDP(7 
NDP(  7 


NDP(1 
NDP(1 
NDP(2 
NDP(2 
NDP(3 
NDP(3 
NDP(4 
NDP(4 
NDP(5 
NDP(5 
NDP(6 


ELSEIFd.EQ.  BIND+1)   THEN 


MESH(  1 , 1)+MESH(  1 , 3)''-DIST 

MESH(1,2)+MESH(1,4)-'^DIST 

MESH(1,1) 

MESH(1,2) 

MESH( PERND, 1)+MESH( PERND, 3 )*DIST 

MESH(  PERND,  2  )+MESH(  PERND,  4)^^DIST 

MESH( PERND,!) 

MESH(PERND,2) 

MESH(PERND+1,1) 

MESH(PERND+1,2) 

MESH(2,1) 

MESH(2,2) 

MESH(2,1)+MESH(2,3)*DIST 

MESH(2,2)+MESH(2,4)''^DIST 


MESH(BIND+3,1)+MESH(BIND+3,3)^'DIST 

MESH(BIND+3,2)+MESH(BIND+3,4)*DIST 

MESH(BIND+3,1) 

MESH(BIND+3,2) 

MESH(PERND+BIND,1) 

MESH(PERND+BIND,2) 

MESH(BIND+1,1) 

MESH(BIND+1,2) 

MESH( BIND+1 , 1)+MESH( BIND+1 , 3 ) -DIST 

MESH(  BIND+1,  2)+MESH(  BIND+1, 4)''^DIST 

MESH(BIND+2,1)+MESH(BIND+2,3)"DIST 
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NDP(6,2)   =  MESH(BIND+2,2)+MESH(BIND+2,4)*DIST 
NDP(7,1)   =  MESH(BIND+2,1) 
NDP(7,2)   =  MESH(BIND+2,2) 
ELSEIFd.EQ.  BIND+2)    THEN 

WRITE(6,-'0  '  ERRORED  OUT  IN  SUBROUTINE  FINDER,  YOU  ATTEMPTED' 
WRITE(6,")   '      TO  CALL  FINDER  WITH  I  =  BIND  +  2! ' 
WRITE(6,*)  •    THIS  ROW  HAS  NO  ELEMENTS  AND  NO  COORDINATES' 
ELSE 
C  NODE  1 

NDP(1,1)  =  MESH(PERND-I+2,1)+MESH(PERND-I+2,3)''^DIST 
NDP(1,2)  =  MESH(PERND-I+2,2)+MESH(PERND-I+2,4)*DIST 
C  NODE  2  TO  THE  BISECTION  SEGMENT 

DO  10,  J  =  2,  N0DES(PERND-I+2) 

NDP(J,1)=MESH(PERND-I+2,1)-(J-2)*(MESH(PERND-I+2,1)- 
C  MESH(PERND+I-l,l))/(N0DES(PERND-I+2)-2) 

NDP(J,2)=MESH(PERND-I+2,2)-(J-2)^KMESH(PERND-I+2,2)- 
C  MESH(PERND+I-l,2))/(N0DES(PERND-I+2)-2) 

10         CONTINUE 

C  BISECTION  SEGMENT  TO  THE  RIGHTS IDE  SURFACE 

DO  20,  J  =  3,  NODES(I) 

NDP(N0DES(PERND-I+2)+J-2,l)=MESH(PERND+I-l,l)+(J-2)''^ 
C  ( MESH( 1,1) -MESH( PERND+I - 1 , 1 ) ) / ( NODES( I ) -2 ) 

NDP(N0DES(PERND-I+2)+J-2,2)=MESH(PERND+I-l,2)+(J-2)''^ 
C  (MESH(I,2)-MESH(PERND+I-l,2))/(N0DES(I)-2) 

20         CONTINUE 
C  Ith  ROWS  LAST  NODE 

NDP(N0DES(PERND-I+2)+N0DES(  I)  -1, 1)=MESH(  1 , 1)+MESH(  1 ,3)'>DIST 
NDP(N0DES(PERND-I+2)+N0DES(I)-l,2)=MESH(I,2)+MESH(I,4)^'-DIST 
C  I+lth  ROWS  FIRST  NODE 

NDP(N0DES(PERND-I+2)+N0DES(I),l)=MESH(PERND-I+l,l)+MESH 
C     (PERND-I  +  1,3)-'^DIST 

NDP(N0DES(PERND-I+2)+N0DES(I),2)=MESH(PERND-I+l,2)+MESH 
C     (PERND-I+1,4)''^DIST 
C  I+ITH  ROW  (NODE  2)  TO  THE  BISECTION  SEGMENT 

DO  30,  J  =  2,  N0DES(PERND-I+1) 

NDP( J+NODES( PERND- 1+2 )+NODES( I ) - 1 , 1 )=MESH( PERND- I+l , 
C  l)-(J-2)-nMESH(PERND-I  +  l,l)-MESH(  PERND+I,  1))/ 

C  (N0DES(PERND-I+l)-2) 

NDP(J+NODES(PERND-I+2)+NODES(I)-l,2)=MESH(PERND-I+l, 
C  2)-(J-2)''^(MESH(PERND-I  +  l,2)-MESH(PERND+I,2))/ 

C  (N0DES(PERND-I+l)-2) 

30         CONTINUE 

C  I+lth  ROW  BISECTION  SEGMENT  TO  THE  RIGHTSIDE  SURFACE 

DO  40,  J  =  3,  N0DES(I+1) 

NDP(  LAND-NODES(  I+l)+J- 1 , 1  )=MESH(  PERND+I ,  1 )+( J-2)'^ 
C  ( MESH( 1+1,1) -MESH( PERND+I , 1 ) ) / ( NODES( I+l ) -2 ) 

NDP(LAND-N0DES(I+l)+J-l,2)=MESH(PERND+I,2)+(J-2)* 
C  (MESH( I+l ,2) -MESH( PERND+I ,2) ) /(NODES ( I+l) -2) 


=  MESH(I+1,1)+MESH(I+1,3)'''DIST 
=  MESH(I+1,2)+MESH(I+1,4)^'DIST 


40 

CONTINUE 

C 

LAST  NODE 

NDP(LAND,1) 

NDP(LAND,2) 

C 

ENDIF 

RETURN 
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10 

c 


END 


SUBROUTINE  VARINT( J, F, ALPHA, BETA, AREA, LND) 

GENERATING  VARIATIONAL  FINITE  ELEMENT  AREA  INTEGRATIONS  OF  THE 
LINEAR  BASIS  FUNCTION  LAGRANGIAN  FOR  THE  HELMHOLTZ  EQUATION. 
THESE  ARE  RETURNED  IN  F(3,3).  X(3)  AND  Y(3)  ARE  THE  WAVENUMBER 
NORMALIZED  CARTESIAN  COORDINATES  OF  THE  TRIANGLE  VERTICES. 
X  =  Ko'-^x,  Y  =  Ko^'^y    ALPHA  AND  BETA  ARE  COMPLEX 
MATERIAL  PARAMETERS  WITHIN  THE  ELEMENT. 

FOR  TM  INCIDENCE:  ALPHA  =  1/ur;  BETA  =  er 
FOR  TE  INCIDENCE:  ALPHA  =  1/er;  BETA  =  ur 

OUTPUTS  :    F(3,3)  -  FINITE  ELEMENT  AREA  INTEGRATION 
AREA   -  AREA  OFF  A  TRIANGLE 

COMPLEX  ALPHA,  BETA,  B12,  F(3,3) 

REAL  NDP(200,2),  X(3),  Y(3),  T(3,3),  AREA,  DET 

INTEGER  L,  K,  J,  LND(0:  200,3) 

COMMON/ BLK5/NDP 

X(l)  =  NDP(LND(J,1),1) 
X(2)  =  NDP(LND(J,2),1) 
X(3)  =  NDP(LND(J,3),1) 
Y(l)  =  NDP(LND(J,1),2) 
Y(2)  =  NDP(LND(J,2),2) 
Y(3)  =  NDP(LND(J,3),2) 
DET  =  X(2)''--Y(3)  +  X(3)- 
CX(1)*Y(3)  -  X(2)^>Y(1) 
AREA  =  ABS(0.5"DET) 
B12  =  BETA/12. 


^Y(l)  +  X(1)^'^Y(2)  -  X(3)''^Y(2) 


T(l,l) 
T(l,2) 
T(l,3) 
T(2,l) 
T(2,2) 
T(2,3) 
T(3,l) 
T(3,2) 
T(3,3) 
DO  10, 
DO 


(Y(2)  -  Y(3))/DET 

(Y(3)  -  Y(1))/DET 

(Y(l)  -  Y(2))/DET 

(X(3)  -  X(2))/DET 

(X(l)  -  X(3))/DET 

(X(2)  -  X(1))/DET 

X(3)^-Y(2))/DET 
X(1)''^Y(3))/DET 
X(2)'-Y(1))/DET 


=  (X(2)^''Y(3) 
=  (X(3)''^Y(1) 
=  (X(l)--'-Y(2) 
K  =  1,  3 
10,  L  =  1,  3 

F(K,L)  =  ALPHA*(T(1,K)*T(1,L)  +  T(2,K)'VT(2,L) )  -  B12 
IF(K.EQ.  L)  F(K,L)  =  F(K,L)  -  B12 
F(K,L)  =  AREA->F(K,L) 

RETURN 
END 


SUBROUTINE  LOADER( BCOND , OFFSET , ALPHA , BETA , NABC , IMX , NBMX , EO , SURBC , 
CCHARl, LINE, MODE, XORIGIN,YORIGIN,SVEC, CHARS) 


COMPLEX  A(50,50),B(50,50),C(50,50),P(50,100) 
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COMPLEX  F(3,3),FROW(100,3,3),BCOND(100),LINE(50) 

COMPLEX  ALPHA,  BETA,  DETERM,  SURBC(IOO),  SVEC(50,100) 

REAL  OFFSET, MINAREA,MAXAREA, AREA, RATIO, E0,K0,XORIGIN,YORIGIN 

REAL  UBCOND 

INTEGER  I,J,JD,K,L,NDT0P,NDB0T,NDT0T,N0D,KND,ND(3),LEL,LAND 

INTEGER  LND(0: 200,3),  NDL(200,4),  NCT(200),  PERND,  BIND,  JJ 

INTEGER  N0DES(200),  MINROW,  MINEL,  MAXROW,  MAXEL,  TCALL,  NL 

INTEGER  N,M,NBMX,NABC(100,3),IN0RM,IMX,M0DE 

CHARACTER'n  CHARl,  CHAR2,  CHAR3 ,  CHAR4,  CHARS 
C 

C0MM0N/BLK2/PERND , BIND 

C0MM0N/BLK3/N0DES 

C0MM0N/BLK4/LND,NDL,NCT 

C0MM0N/BLK7/MINAREA, MINROW, MINEL, MAXAREA, MAXROW, MAXEL, AREA 

C0MM0N/BLK8/A,B,C,P 

C0MM0N/BLK9/CHAR2,  CHAR3 ,  CHAR4 
C 

UBCOND  =1.0 

I NORM  =  0 

IMX  =  BIND  +  2 

I  =  1 

TCALL  =  BIND  +  1 
C 

IF((CHAR3.EQ.  'D').OR.  (CHAR3.EQ.  'd'  ))  THEN 
WRITE(20,1030) 
WRITE(20,1040) 
WRITE(20,1050) 
WRITE(20,1060) 
WRITE(20,1070) 

END  IF 

WRITE(''S1000)  I, TCALL 
C 

CALL  SORTER(I,LEL,LAND) 

CALL  FINDERC I, LAND, OFFSET) 

IF(.NOT.  (( CHARS.  EQ.  'M').0R.  (CHARS.  EQ.  'm')))  THEN 

CALL  BNDC( I, BCOND,E0,SURBC, CHARl, ALPHA, BETA, LINE, MODE, XORIGIN, 
CY0RIGIN,CHAR4) 

END  IF 

CALL  FILL( I, LEL,FROW, ALPHA, BETA) 

IF(.NOT.  ((CHARS.  EQ.  'M').0R.  (CHARS.  EQ.  'm'  )))  THEN 

CALL  ZERO 

END  IF 
C     ESTABLISH  THE  NUMBER  OF  NODES  IN  THE  I  =  1  AND  I  =  2  ROWS 

NDTOP  =  2 

NDBOT  =  S 

NDTOT  =  7 
C 

C  B,C&PF0RI=1    (TOP) 

C 

J  =    1 

DO   70,    JD  =   1,    NCT(2) 
L  =  NDL(2,JD) 
DO  50,    K  =   1,    3 

ND(K)  =  LND(L,K) 
IF(ND(K).EQ.  2)  KND  =  K 
50         CONTINUE 
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DO  60,  K  =  1,  3 
NL  =  ND(K) 
C     BOUNDARY  CONDITION 

IF(NL.  EQ.  1)  THEN 

P(J,1)  =  -UBCOND''^FROW(L,KND,K)  +  P(J,1) 
C     UPPER  ROW 

ELSEIF(NL.  EQ.  2)  THEN 

B(J,1)  =  FROW(L,KND,K)  +  B(J,1) 
C     LOWER  ROW 

ELSE IF(  ( NL.  GT.  3 ) .  AND.  ( NL.  LT.  7 )  )  THEN 

C(J,NL-3)  =  FROW(L,KND,K)  +  C(J,NL-3) 
C     ERROR 

ELSE 

WRITEC*,*)  NL, 'ERROR  -  INDEX  EQUALS  3  OR  7' 
ENDIF 
60        CONTINUE 
70    CONTINUE 
C 

IF((CHAR2.EQ.  'l').OR.  (CHAR2.EQ.  'i'))  THEN 
WRITEd,*)  I 
DO  72,  M  =  1,  NABC(I,2) 

WRITE(1,1020)  (REAL(B(M,N)),  N  =  1,  NABC(I,2)) 

72  CONTINUE 

WRITECl,''^)  '  ' 
DO  73,  M  =  1,  NABC(I,2) 

WRITE(1,1020)  (REAL(C(M,N)),  N  =  1,  NABC(I,3)) 

73  CONTINUE 

WRITECl,'-)  '  ' 
DO  74,  M  =  1,  NABC(I,2) 

WRITE(1,1020)  (REAL(P(M,N)),  N  =  1,  PERND) 


74 

CONTINUE 

C 

WRITE(1,''0    '    ' 

WRITEd,^'-)    '     ' 

WRITE(1,''0    '    ' 

ENDIF 
C 

IF(.NOT.  (( CHARS.  EQ.  'M').0R.  (CHARS.  EQ.  'm')))  THEN 
CALL  MARCH( I , IMX,NBMX,NABC( I , 1) ,NABC( 1,2) ,NABC( 1,3), INORM,SVEC) 
CALL  ZERO 
ENDIF 
C     BEGIN  LOOPING 

DO  900,  1=1,  BIND 
C 

DO  400,  J  =  1,  NDBOT-2 

NOD  =  J  +  NDTOP  +  1 
DO  300,  JD  =  1,  NCT(NOD) 
L  =  NDL(NOD,JD) 
DO  100,  K  =  1,  3 

ND(K)  =  LND(L,K) 
IF(ND(K).EQ.  NOD)  KND  =  K 
100  CONTINUE 

DO  200,  K  =  1,  3 
NL  =  ND(K) 
C     BOUNDARY  CONDITION 

IF((I.EQ.  1).AND.  (NL.  EQ.  D)  THEN 
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P(J,NL)  =  -UBCOND''^FROW(L,KND,K)  +  P(J,NL) 
ELSEIF(NL.  EQ.  1)  THEN 

P(J,PERND-I+2)=-UBC0ND-''^FR0W(L,KND,K)+P(J,PERND-I+2) 
C     SPECIAL  CASE  FOR  NODE  //2  AND  I  =  1 
ELSEIF((I.EQ.  1).  AND.  (NL.  EQ.  2))THEN 

A(J,1)  =  FROW(L,KND,K)  +  A(J,1) 
ELSEIF(NL.  EQ.  NDTOP)  THEN 

P(J,I)  =  -UBCOND^--FROW(L,KND,K)  +  P(J,I) 
ELSEIF(NL.  EQ.  NDTOP+1)  THEN 

P(J,PERND-I+1)  =  -UBC0ND^'-FR0W(L,KND,K)+P(J,PERND-I+1) 
ELSEIF(NL.  EQ.  NDTOT)  THEN 

P(J,I+1)  =  -UBCOND^^FROW(L,KND,K)  +  P(J,I+1) 
C     UPPER  ROW 

ELSEIF(NL.LT.  NDTOP)  THEN 

A(J,NL-1)  =  FROW(L,KND,K)  +  A(J,NL-1) 
C     LOWER  ROW 

ELSEIF(NL.LT.  NDTOT)  THEN 

B(J,NL-NDT0P-1)  =  FROW(L,KND,K)  +  B( J,NL-NDT0P-1) 
C     ERROR 
ELSE 

WRITE('^''0  NL, 'ERROR  -  DUE  TO  INDEX  GREATER  THAN  NODE  TOTAL' 
END  IF 
C 

200  CONTINUE 

300  CONTINUE 

400  CONTINUE 

C 

C     INDEX  FOR  NEXT  ROW  AND  COMPLETE   B,  C,  P  FILLS 
C 

JJ  =  I  +  1 

WRITE(6,1010)  JJ,TCALL 
CALL  SORTER(JJ,LEL,LAND) 
CALL  FINDER(JJ, LAND, OFFSET) 

IF(.NOT.  ((CHAR5.EQ.  'M').0R.  (CHARS.  EQ.  'm')))  THEN 
CALL  BNDC(JJ,BC0ND,E0,SURBC,CHAR1, ALPHA, BETA, LINE, MODE, 
C     X0RIGIN,Y0RIGIN,CHAR4) 
ENDIF 

CALL  FILL( JJ , LEL , FROW , ALPHA , BETA) 
C     ESTABLISH  THE  NUMBER  OF  NODES  IN  THE  I+lth  AND  I+2th  ROWS 
IF(JJ. NE. (BIND+1))  THEN 

NDTOP  =  N0DES(PERND+2-JJ)  +  NODES(JJ)  -  1 
NDBOT  =  N0DES(PERNT)+1-JJ)  +  N0DES(JJ+1)  -  1 
NDTOT  =  NDBOT  +  NDTOP 
ELSE 

NDTOP  =  5 
NDBOT  =  2 
NDTOT  =  7 
ENDIF 
C 

DO   800,    J  =   1,    NDTOP-2 
NOD  =  J  +   1 

DO   700,    JD  =   1,    NCT(NOD) 
L  =  NDL(NOD,JD) 
DO  500,    K  =   1,    3 

ND(K)   =  LND(L,K) 
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IF(ND(K).EQ.  NOD)  KND  =  K 
500  CONTINUE 

DO  600,  K  =  1,  3 
NL  =  ND(K) 
C     BOUNDARY  CONDITION 
IF(NL.  EQ.  1)  THEN 

P(J,PERND-JJ+2)=-UBC0ND*FR0W(L,KND,K)+P(J,PERND-JJ+2) 
ELSEIFCNL. EQ.  NDTOP)  THEN 

P(J,JJ)  =  -UBCOND*FROW(L,KND,K)  +  P(J,JJ) 
ELSEIF(NL.  EQ.  NDTOP+1)  THEN 

P(  J ,  PERND- JJ+1  )  =  -UBCOND''^FROW(  L,KND  ,K)+P(  J ,  PERND- JJ+1) 
ELSEIF((JJ.  EQ.  (BIND+1)),  AND.  (NL.  EQ.  NDTOT)  )  THEN 

C(J,1)  =  FROW(L,KND,K)  +  C(J,1) 
ELSEIFCNL.  EQ.  NDTOT)  THEN 

P(J,JJ+1)  =  -UBCOND''^FROW(L,KND,K)  +  P(J,JJ+1) 
C     UPPER  ROW 

ELSEIFCNL.  LT. NDTOP)  THEN 

BCJ,NL-1)  =  FROWCL,KND,K)  +  BCJ,NL-1) 
C     LOWER  ROW 

ELSEIFCNL. LT.  NDTOT)  THEN 

CCJ,NL-NDT0P-1)  =  FROWCL,KND,K)  +  CC J,NL-NDT0P-1) 
C     ERROR 
ELSE 

WRITER'S ''^)  NL, 'ERROR  -  DUE  TO  INDEX  GREATER  THAN  NODE  TOTAL' 
ENDIF 


C 

600 

CONTINUE 

700 

CONTINUE 

800 

CONTINUE 

C 

IFCCCHAR2.EQ.  'l').OR.  CCHAR2.EQ.  'i'))  THEN 
WRITECl,")  JJ 
DO  91,  M  =  1,  NABCCJJ,2) 

WRITEC1,1020)  CREALCACM,N)),  N  =  1,  NABCCJJ,1)) 

91  CONTINUE 

WRITECl,")  '  ' 
DO  92,  M  =  1,  NABCCJJ,2) 

WRITEC1,1020)  CREALCBCM,N)),  N  =  1,  NABC(JJ,2)) 

92  CONTINUE 

WRITECl,")  '  ' 
DO  93,  M  =  1,  NABCCJJ,2) 

WRITEC1,1020)  CREALCCCM,N)),  N  =  1,  NABCCJJ,3)) 

93  CONTINUE 

WRITECl,^--)  '  ' 
DO  94,  M  =  1,  NABCCJJ,2) 

WRITEC1,1020)  CREALCPCM,N)),  N  =  1,  PERND) 


94 

CONTINUE 

C 

WRITECl,*)  '  ' 

WRITECl,*)  '  ' 

WRITECl,'"^)  '  ' 

ENDIF 

IFC.NOT.  CCCHAR5.EQ.  'M').0R.  CCHAR5.EQ.  'm')))  THEN 
CALL  MARCHCJJ,IMX,NBMX,NABCCJJ,1),NABCCJJ,2),NABCCJJ,3),IN0RM, 
CSVEC) 
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CALL  ZERO 
END  IF 
C 

900   CONTINUE 
C 

C     LOAD  A,  B  &  P  FOR  I  =  BIND  +  2  (BOTTOM) 
C 

J  =  1 

DO  30,  JD  =  1,  NCT(7) 
L  =  NDL(7,JD) 
DO  10,  K  =  1,  3 

ND(K)  =  LND(L,K) 
IF(ND(K).EQ.  7)  KND  =  K 
10  CONTINUE 

DO  20,  K  =  1,  3 
NL  =  ND(K) 
C     BOUNDARY  CONDITION 
IF(NL.  EQ.  6)  THEN 

P(J,JJ+1)  =  -UBCOND-'^FROW(L,KND,K)  +  P(J,JJ+1) 
C     UPPER  ROW 

ELSEIF((NL.  GT.  1).  AND.  (NL.  LT.  5))  THEN 

A(J,NL-1)  =  FROW(L,KND,K)  +  A(J,NL-1) 
C     LOVER  ROW 

ELSEIF(NL.  EQ.  7)  THEN 

B(J,1)  =  FROW(L,KND,K)  +  B(J,1) 
C     ERROR 
ELSE 

WRITE('S'-)  NL, 'ERROR  -  INDEX  EQUAL  TO  1  OR  5* 
ENDIF 
C 

20         CONTINUE 
30    CONTINUE 
C 

IF((CHAR2.EQ.  •l').OR.  (CHAR2.EQ.  'i'))  THEN 
WRITE(1,''^)  TCALL+1 
DO  61,  M  =  1,  NABC(TCALL+1,2) 

WRITE( 1,1020)  (REAL(A(M,N)),  N  =  1,  NABC ( TCALL+1 , 1) ) 

61  CONTINUE 

WRITE(1,'>)  '  ' 
DO  62,  M  =  1,  NABC (TCALL+1, 2) 

WRITE(1,1020)  (REAL(B(M,N)),  N  =  1,  NABC(TCALL+1 ,2) ) 

62  CONTINUE 

WRITE(1,''0  '  ' 
DO  64,  M  =  1,  NABC (TCALL+1, 2) 

WRITE( 1,1020)  (REAL(P(M,N)),  N  =  1,  PERND) 
64         CONTINUE 
C 

WRITEd,''-)  '  ' 
WRITEd,'"^)  '  ' 
WRITE(1,^>-)  '  ' 
ENDIF 
C 

WRITE(6,''0  '     FINAL  INVERSION' 
IF(.NOT.  ( (CHARS.  EQ.  'M').0R.  (CHARS.  EQ.  'm')))  THEN 
CALL  MARCH(JJ+1,IMX,NBMX,NABC(JJ+1,1),NABC(JJ+1,2),NABC(JJ+1,3), 
CINORM,SVEC) 
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END  IF 


C 

1000 
1010 
1020 
1030 
1040 
1050 
1060 
1070 
1080 
1090 
1100 

c 


END 


MINIMUM  AREA  =  ' ,MINAREA 

AT  ROW  ',MINROW  '     AND  ELEMENT  NUMBER  ' ,MINEL 

MAXIMUM  AREA  =   jMAXAREA 

AT  ROW  ',MAXROW,'      AND  ELEMENT  NUMBER  ' ,MAXEL 


WRITE ( 2, ''0  '  END  END 

WRITE  ( 6, --'0 

WRITE  ( 6, ''^) 

WRITE(6,'^) 

WRITE(6,^'0 

WRITE(6,''0 

RATIO  =  MAXAREA/MINAREA 

WRITE(6,-'"^)  'AREA  RATIO   =  ', RATIO 

IF( RATIO.  GT.  2.5)  THEN 

WRITE (6,*)  'YOU  SHOULD  CONSIDER  ABORTING  THIS  RUN  AND 
WRITE(6,''0  'LOOKING  AT  THE  MESH  IN  CURVE  DIGITIZER  ' 
WRITECe,'^)  'A  BETTER  METHOD  MAY  BE  AVAILABLE  ' 

ENDIF 

IF( ( CHARS.  EQ.  'D').0R.  (CHAR3.EQ.  'd'))  THEN 
WRITE(20,1080) 
WRITE(20,1090) 
WRITE(20,1100) 

ENDIF 


FORMATC 
FORMAT( 
FORMATC 
FORMATC 
FORMATC 
FORMATC 
FORMATC 
FORMATC 
FORMATC 
FORMATC 
FORMATC 

RETURN 
END 


16,'  OUT  OF' ,16,'  CALLS  ') 
16,'  OUT  OF' ,16,'  CALLS  ') 
20CE14.  8,1X,E14.  8, IX)) 
6X,'CALL  COMPRS' ) 
6X,'CALL  NOBRDR') 
6X,'CALL  PAGEC8.  0,10.  0)') 
6X,'CALL  AREA2DC5.  0,7.  0)') 
6X,'CALL  FRAME') 
6X,'CALL  DONEPL') 
6X,'ST0P') 
6X,'END') 


SUBROUTINE  FILLC I , LEL , FROW , ALPHA , BETA) 

COMPLEX  FC3,3),  FROWC 100 , 3 , 3) ,  ALPHA,  BETA,  A,  B 

REAL  AREA,  MINAREA,  MAXAREA 

REAL  NDPC200,2) 

INTEGER  I,  J,  K,  L,  LEL,  LNDCO: 200 , 3) ,  NDLC200,4),  NCT(200) 

INTEGER  MINROW,  MINEL,  MAXROW,  MAXEL,  M 

CHARACTER" 1  CHAR2 ,  CHAR3 ,  CHAR4 

C0MM0N/BLK4/LND,NDL,NCT 

C0MM0N/BLK5/NDP 

C0MM0N/BLK7/MINAREA, MINROW, MINEL, MAXAREA, MAXROW, MAXEL, AREA 

C0MM0N/BLK9/CHAR2,  CHAR3,  CHAR4 

IFCCCHAR2.EQ.  'l').OR.  CCHAR2.EQ.  'i'))  THEN 

WRITECll,'"')  'ROW  NUMBER  =  ',1 
ENDIF 
DO  30,  J  =  1,  LEL 

IFC  C  CHAR4.  EQ.  '  U '  ) .  OR.  C  CHAR4.  EQ.  '  u '  )  )  THEN 
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A  =  ALPHA 
B  =  BETA 
ELSE 

IF((J.  EQ.  1).0R.  (J.EQ.  2).0R.  ( J.  EQ.  LEL-1).  OR.  (J.EQ.  LEL)) 
C  THEN 

A  =  (1.0,0.0) 
B  =  (1.  0,0.  0) 
ELSE 

A  =  ALPHA 
B  =  BETA 
ENDIF 
ENDIF 

CALL  VARINT(J,F,A,B,AREA,LND) 
IF((J.  GT.  2).  AND.  ( J.  LT.  LEL-1) )  THEN 
IF(AREA. LT. MINAREA)  THEN 
MINAREA  =  AREA 
MINROW  =  I 
MINEL  =  J 
ELSEIF( AREA. GT.  MAXAREA)  THEN 
MAXAREA  =  AREA 
MAXROW  =  I 
MAXEL  =  J 
ENDIF 
ENDIF 
DO  20,  K  =  1,  3 

WRITE(2,^0  NDP(LND(J,K),1),  NDP(LND( J,K) ,2) 
DO  10,  L  =  1,  3 

FROW(J,K,L)  =  F(K,L) 
10  CONTINUE 

IF((CHAR2.EQ.  'r).OR.  (CHAR2.EQ.  'i'))  THEN 

WRITE(11,1000)  J,K,(REAL(F(K,L)),  L  =  1,3) 
ENDIF 
20         CONTINUE 

IF(  ( CHAR2.  EQ.  '  I '  ) .  OR.  ( CHAR2.  EQ.  '  i '  ) )  THEN 

WRITE ( 11, -'0  '  ' 
ENDIF 

WRITE(2,''0  NDP(LND(J,1),1),  NDP(LND(  J,  1)  ,2) 
WRITE(2,")  '  999990   999990  ' 
C     DISSPLA  PROGRAM  GENERATION 

IF((CHAR3.EQ.  'D').0R.  (CHAR3.EQ.  'd'))  THEN 

WRITE(20,1010)  NDP(LND(J,1),1),  NDP( LND( J, 1) , 2) 
WRITE(20,1020)  NDP(LND( J,2) , 1) ,  NDP(LND( J,2) ,2) 
WRITE(20,1020)  NDP(LND( J,3) , 1) ,  NDP(LND( J,3) ,2) 
WRITE(20,1020)  NDP(LND( J, 1) , 1) ,  NDP(LND( J, 1) ,2) 
ENDIF 
C 

30    CONTINUE 
C 

1000  F0RMAT(1X,2(I3,2X),3X,3(F8.  5,2X)) 
1010  F0RMAT(6X,'CALL  STRTPT( ' ,F8. 5 , ' , ' ,F8. 5 , ' ) ' ) 
1020  F0RMAT(6X,'CALL  C0NNPT( ' ,F8. 5 , ' , ' ,F8.  5, ' ) ' ) 
C 

RETURN 
END 
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SUBROUTINE  ZERO 

ZERO  A,  B,  C,  AND  P  MATRICES 

COMPLEX  A(50,50),  B(50,50),  C(50,50),  P(50,100) 

INTEGER  J,  K,  L 

C0MM0N/BLK8/A,B,C,P 


DO  50,  J  =  1,  50 

DO  40,  K  =  1,  50 

A(J,K)  =  CMPLX(0.  0,0.0) 
B(J,K)  =  CMPLX(0.  0,0.  0) 
C(J,K)  =  CMPLX(0.  0,0.  0) 
40         CONTINUE 

DO  60,  L  =  1,  100 

P(J,L)  =  CMPLX(0.0,0.  0) 
60         CONTINUE 

50    CONTINUE 

RETURN 

END 
C 
C 

SUBROUTINE  BNDC( I , BCOND , EO , SURBC , CHARl , ALPHA , BETA , LINE , MODE , 
CXORIGIN , YORIGIN , CHAR4) 
C 

C     BOUNDARY  CONDITION  FILL  FOR  A  PLANE  WAVE 
C     EO  =  FIELD  STRENGTH 
C     KO  =  WAVE  NUMBER  (2*PI/WAVELENGTH) 
C 

C     BOUNDARY  CONDITION  FILL  FOR  A  CYLINDRICAL  BOUNDARY  CONDITION 
C     EO  =  FIELD  STRENGTH 

C     MODE  =  MODE  NUMBER  FOR  CYLINDRICAL  BOUNDARY  CONDITIONS 
C 

COMPLEX  SURBC(IOO),  VALUE(50),  BCOND(IOO),  ALPHA,  BETA,  LINE(50) 

COMPLEX  K,  RAK 

REAL  NDP(200,2),  EO,  KR,  KI ,  XORIGIN,  YORIGIN,  RA,  RB 

INTEGER  I,  J,  NDTOP,  NDBOT,  NDTOT,  NODES(200),  PERND,  BIND,  MODE 

CHARACTER" 1  CHARl,  CHAR4 
C 

COMMON/ BLK2/ PERND , BIND 

C0MM0N/BLK3/N0DES 

C0MM0N/BLK5/NDP 


IF((CHAR1.EQ.  'P').0R.  (CHARl.  EQ.  p'))  THEN 

IF(  (  CHAR4.  EQ.  '  U '  ) .  OR.  (  CHAR4.  EQ.  '  u '  )  )  THEN 
KR  =  REAL(CSQRT(BETA)) 
KI  =  ABS(AIMAG(CSQRT(BETA))) 
ELSE 

KR  =  1.  0 
KI  =  0.0 
END  IF 
C 

C     CHECK  ON  WHETHER  WE  ARE  USING  ER  =  ALPHA  OR  BETA 
C 

IF(I.EQ.  1)  THEN 

BCOND(l)  =  E0-'^EXP(KI^^NDP(1,2))''>-CMPLX(COS(KR''^NDP(1,2)), 
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C     SIN(KR''^NDP(1,2))) 

BCOND(  PERND)  =  EO*EXP(KI'>NDP(  3 , 2)  )*CMPLX(COS(KR*NDP(  3,2)), 
C     SIN(KR-"-NDP(3,2))) 

BC0ND(2)  =  EO'>EXP(KI'VNDP(7,2))*CMPLX(COS(KR->NDP(7,2)), 
C     SIN(KR-"-NDP(7,2))) 

VALUE(l)  =  E0--'^EXP(KI'^NDP(2,2))''-CMPLX(COS(KR*NDP(2,2)), 
C     SIN(KR''-NDP(2,2))) 

WRITE(10,--^)  BCOND(l) 

WRITE(10,1000)  VALUE(l) 

WRITEdO,*)  '  ' 

SURBC(l)  =  VALUE(l) 
ELSEIF(I.LE.BIND)  THEN 

NDTOP  =  NODES(I)  +  N0DES(PERND-I+2)  -  1 

NDBOT  =  N0DES(I+1)  +  NODES( PERND -I+l)  -  1 

NDTOT  =  NDTOP  +  NDBOT 

BC0ND(PERND-I+1)  =  EO'^^EXPCKI^NDPCNDTOP+l ,  2)  )*CMPLX(COS(KR* 
C     NDP(  NDTOP+1 ,  2  )  )  ,  SIN(  KR''-NDP( NDTOP+1 ,  2 )  )  ) 

BC0ND(I+1)  =  E0-''fEXP(KI''^NDP(NDTOT,2))*CMPLX(COS(KR*NDP(NDTOT, 
C     2)),SIN(KR'^NDP(NDT0T,2))) 

WRITE(  10, -'0  BC0ND(PERND-I+2) 

DO  10,  J  =  2,  NDTOP- 1 

VALUE(J)  =  E0-"EXP(KI''^NDP(J,2))*CMPLX(COS(KR''fNDP(J,2)), 
C  SIN(KR-'^NDP(J,2))) 

IF(J.  EQ.  2)  THEN 

SURBC( PERND -1+2)  =  VALUE(2) 
ELSEIF(J.EQ.  (NDTOP- 1))  THEN 

SURBC(I)  =  VALUE( NDTOP- 1) 
ENDIF 
10         CONTINUE 

LINE(I)  =  VALUEC (NDTOP+1 )/2) 

WRITE( 10,1000)  (VALUE(J),  J  =  2,  NDTOP- 1) 

WRITEdO,")  BCOND(I) 

WRITE(10,''O  '  ' 
C 

ELSEIF(I.EQ.  BIND+1)  THEN 

BC0ND(I+1)  =  E0--'-EXP(KI^^NDP(6,2))"CMPLX(C0S(KR'-NDP(6,2)), 
C     SIN(KR-'^NDP(6,2))) 

VALUE(2)  =  E0^--EXP(KI*NDP(2,2))^^CMPLX(C0S(KR"NDP(2,2)), 
C     SIN(KR''^NDP(2,2))) 

VALUE(3)  =  E0-''-EXP(KI*NDP(3,2))^fCMPLX(COS(KR"NDP(3,2)), 
C     SIN(KR--'-NDP(3,2))) 

VALUE(4)  =  E0"EXP(KI'''NDP(4,2))^-CMPLX(C0S(KR-^NDP(4,2)), 
C     SIN(KR-'>NDP(4,2))) 

WRITE( 10,^0  BC0ND(I+2) 

WRITEdO,  1000)  VALUE(2),  VALUE(3),  VALUE(4) 

SURBC(BIND+3)  =  VALUE(2) 

SURBC( BIND+1)  =  VALUE(4) 

LINE(I)  =  VALUE(3) 

WRITEdO,*)  BCOND(I) 

WRITEdO,''-)  '  ' 

VALUE(2)  =  E0'VEXP(KI'mDP(7,2))*CMPLX(C0S(KR*NDP(7,2)), 
C     SIN(KR--NDP(7,2))) 

WRITEdO,  1000)   VALUE(2) 

WRITEdO, -'0    BC0ND(I+1) 

SURBC(BIND+2)   =  VALUE(2) 
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ENDIF 
C 

ELSEIF((CHAR1.EQ.  'C').OR.  (CHARl.EQ.  'c'))  THEN 
C 

IF(I.EQ.  1)  THEN 

RA  =  SQRT((NDP(2,1)-X0RIGIN)''"V2+(NDP(2,2)-Y0RIGIN)**2) 
RB  =  SQRT((NDP(1,1)-X0RIGIN)'V-V2+(NDP(1,2)-Y0RIGIN)''"^2) 
WRITE (--'s^'O  BETA 
K  =  CSQRT(BETA) 
RAK  =  RA-'^K 

CALL  CYLBC(RA,RB,RAK,NDP(l,l),NDP(l,2),XORIGIN,YORIGIN,E0, 
C     BC0ND(1),SURBC(1),K) 

CALL  CYLBC(RA,RB,RAK,NDP(3,l),NDP(3,2),XORIGIN,YORIGIN,E0, 
C     BCOND(PERND),SURBC(PERND),K) 

CALL  CYLBC(RA,RB,RAK,NDP(7,l),NDP(7,2),XORIGIN,YORIGIN,EO, 
C     BCOND(2),SURBC(2),K) 

WRITE(10,''O  BCOND(l),  SURBC(l) 
WRITE(10,''O  '  ' 
ELSEIF(I.LE.BIND)  THEN 

NDTOP  =  NODES(I)  +  N0DES(PERND-I+2)  -  1 
NDBOT  =  N0DES(I+1)  +  NODES( PERND-I+1)  -  1 
NDTOT  =  NDTOP  +  NDBOT 

CALL  CYLBC(RA,RB,RAK,NDP(NDT0P+1,1),NDP(NDT0P+1,2),X0RIGIN, 
C     YORIGIN,EO,BCOND(PERND-I+1),SURBC(PERND-I+1),K) 

CALL  CYLBC(RA,RB,RAK,NDP(NDT0T,1),NDP(NDT0T,2),X0RIGIN, 
C     Y0RIGIN,E0,BC0ND(I+1),SURBC(I+1),K) 

WRITE(10,''O    BC0ND(PERND-I+2),    SURBC(PERND-I+2) 
WRITECIO,^'-)    BCOND(I),    SURBC(I) 
WRITE(10,'>)    '     ' 
ELSEIFd.EQ.  BIND+1)   THEN 

WRITEdO,''-)    BC0ND(I+2),    SURBC(I+2) 
WRITEdO,^-)    BCOND(I),    SURBC(I) 
WRITEdO,''-)    '     ' 

CALL  CYLBC(RA,RB,RAK,NDP(6,1),NDP(6,2),X0RIGIN, 
C     YORIGIN,EO,BCOND(I+1),SURBC(I+1),K) 
WRITE(10,''O  BCONDd+1),  SURBC(I+1) 
WRITE(10,-'O  '  ' 
ENDIF 

^  N 

ELSE         ) 
RETURN 

ENDIF 
C 

1000  FORMAT(1X,50(E14.  8,2X,E14.  8,2X)) 
C 


RETURN 
END 


SUBROUTINE  CYLBC(RA,RB,RAK,X,Y,XORIGIN,YORIGIN,E0,BC,PSI ,K) 

COMPLEX  SQRTMl,  JORB,  JIRB,  JORAK,  J IRAK,  HORA,  HIRA,  HORB,  HIRB 
COMPLEX  DELTAN,  AN,  PSI,  JORA,  JIRA,  K,  RAK,  BC 
REAL  PI,  XORIGIN,  YORIGIN,  X,  Y,  RA,  RB ,  EO,  PHI 

PI  =  4.  0^--ATAN(l.  0) 
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SQRTMl  =  CMPLX(0.  0,1.  0) 

PHI  =  ATAN2(X  -  XORIGIN,  Y  -  YORIGIN) 

CALL  BES1(CMPLX(RA,0. 0) , J0RA,J1RA) 

CALL  BES1(CMPLX(RB,0. 0) , JORB , JIRB) 

CALL  BES1(RAK,J0RAK,J1RAK) 

CALL  HAN1(CMPLX(RA,0. 0) ,H0RA,H1RA) 

CALL  HAN1(CMPLX(RB,0.0),H0RB,H1RB) 
C 

DELTAN  =  J1RB-V(J1RAK-''^(H0RA-H1RA/RA)-K*(J0RAK-J1RAK/RAK)'^H1RA) 
CH1RB*(J1RAK-''KJ0RA-J1RA/RA)-K-'HJ0RAK-J1RAK/RAK)*J1RA) 

AN  =  -2.0-'^SQRTMl/(PI''^RA''^DELTAN) 

BC  =  EO''^COS(PHI) 

PSI  =  AN'>JlRAK-''-BC 
C 

RETURN 

END 
C 

c 

SUBROUTINE  BES1( Z , JO, Jl) 
C 

C       Computing  Bessel  Functions  for  n  =  0,  1  with 
C       Complex  Argument  Z.    Direct  Power  Series  Method  for 
C       CABS(Z)  . LE.  6   and  Hankel's  Asymptotic  Formula  for 
C       CABS(Z)  .GT.  6. 
C       Written  11/5/87  by  M. A.  Morgan 
C 

INTEGER  M,M2 

REAL  C(34),DM,F(34),G0,P(34),Pi,P2 

COMPLEX  Z,Z2,Z3,Z4,J0,J1,AM,CL,P0,P1,Q0,Q1,C0,C1,S0,S1 
Pi=3. 1415927 
P2=2.  0/PI 

IF(CABS(Z).LE.6.0)  THEN 
C 

C   Utilizing  the  Direct  Power  Series  Method 
C 

G0=  1. 781072 
Z2=0.  5-'-Z 
CL=CLOG(G0^^Z2) 
C 

C       Computing  F(m)  =  m  !   and  P(m)  =  1  +  1/2  +  1/3  +  +  1/m 

C 

F(l)  =  1.0 
P(l)=l.  0 
DO  11  M=2,34 

F(M)=M'''F(M-1) 
P(M)=P(M-1)+1.0/M 
11  CONTINUE 

C 

C       Computing  Power  Series  Coefficients 
C 

DM=-1.0 

DO  22  M=l,34 

C(M)=DM/(F(M)'>F(M)) 
DM=-DM 
22  CONTINUE 

C 
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33 


Computing  JO  and  Jl 

J0=(1.  ,0.  ) 
J1=(0.  ,0.  ) 
M=0 
M=M+1 

AM=C(M)^KZ2**M2) 
JO=JO+AM 
J1=J1-M-"AM 

IF((CABS(AM).GT.  1.0E-10).AND.  (M.  LT.  34))  GO  TO  33 
J1=J1/Z2 
return 
ELSE 

Hankel'  Asymptotic  Formula  (Abram.  &  Stegun  p.  364) 

Z2=Z^''Z 

Z3=Z^'^Z2 

Z4=Z^'>-Z3 

P0=1. 0-. 0703125/Z2+.  1121521/Z4 

Q0=-. 125/Z+.  0732422/Z3 

Pl=l.  0+.  1171875/Z2-.  1441956/Z4 

Ql=. 375/Z-.  10253906/Z3 

CO=CCOS(Z-.  25'>PI) 

S0=CSIN(Z-.25"PI) 

C1=CC0S(Z-.  75-'-PI) 

S1=CSIN(Z-.  75--PI) 

AM=CSQRT(P2/Z) 

JO=AM"(PO--'--CO-QO''-SO) 

Jl=AM-"-(Pl''^Cl-Ql''^Sl) 
END  IF 
RETURN 
END 


SUBROUTINE  HAN1(Z,H0,H1) 

Computing  Hankel  Functions  for  n  =  0,  1  with 
Complex  Argument,  Z.    Direct  Power  Series  Method  for 
CABS(Z)  . LE.  5   and  Hankel 's  Asymptotic  Formula  for 
CABS(Z)  .GT.  5.    Written  11/6/87  by  M. A.  Morgan 

INTEGER  M,M2 

REAL  C(34),DM,F(34),G0,P(34),Pi,P2 

COMPLEX  Z,Z2,Z3,Z4,J0,J1,Y0,Y1,AM,CL,P0,P1,Q0,Q1 

COMPLEX  EO,El,XO,Xl,HO,Hl,j 

PI=3.  1415927 

P2=2.0/PI 

j=(0.  ,1.) 

IF(CABS(Z).LE.5.0)  THEN 

Direct  Power  Series  Method 

G0=  1.78072 
Z2=0.  5"Z 
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CL=CL0G(G0-'-Z2) 
C 

C   Computing  F(m)  =  rn  !   and  P(ra)  =  1  +  1/2  +  1/3  + +  1/m 

C 

F(l)  =  l.  0 

P(l)  =  1.0 

DO  11  M=2,34 

F(M)=MV^F(M-1) 
P(M)=P(M-1)  +  1.  0/M 
11  CONTINUE 

C 

C   Computing  Power  Series  Coefficients 
C 

DM=-1.  0 

DO  22  M=l,34 

C(M)=DM/(F(M)*F(M)) 
DM=-DM 
22  CONTINUE 

C 

C   Computing  JO  and  Jl 
C 

J0=(1.  ,0.  ) 

J1=(0.  ,0.  ) 

M=0 
33  M=M+1 

M2=2-->M 

AM=C(M)'nZ2^'">M2) 

JO=JO+AM 

Jl=Jl-M-"-AM 

IF((CABS(AM).GT.  1.0E-10).AND.  (M.  LT.  34))  GO  TO  33 

J1=J1/Z2 
C 

C   Computing  YO  and  Yl 
C 

M=0 

YO=CL"JO 

Yl=Z2''^CL''ai-0.5'''J0 
44  M=M+1 

M2=2-'^M 

AM=C  (  M  )  ^--P (  M  )  ''K  Z 2^>-''-M2  ) 

YO=YO-AM 

Yl=Yl+M--'-AM 

IF((CABS(AM).GT.  l.OE-10).  AND.  (M.  LT.  34))  GO  TO  44 

Y0=P2"Y0 

Y1=P2--VY1/Z2 

HO=JO-j^'-YO 

Hl=Jl-j''^Yl 

RETURN 
ELSE 
C 

C   Hankel'  Asymptotic  Formula  (Abram.  &  Stegun  p.  364 
C 

Z2=Z'>Z 

Z3=Z--'-Z2 

Z4=Z-"-Z3 

P0=1. 0-. 0703125/Z2+. 1121521/Z4 
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Q0=-. 125/Z+.  0732422/Z3 
Pl=l.  0+.  1171875/Z2-.  1441956/Z4 
Ql=.  375/Z-.  10253906/Z3 
X0=(Z-.25*PI) 
X1=(Z-.  75-'^PI) 
EO=CEXP(-j'VXO) 
El=CEXP(-j'>Xl) 
AM=CSQRT(P2/Z) 
HO=AM-'HPO-j''^QO)-''^EO 
Hl=AM-'-(Pl-j^'-Ql)->El 
ENDIF 
RETURN 
END 

SUBROUTINE  MARCH( I , IMX,NBMX,NA,NB,NC, INORM,SVEC) 
C     THIS  ROUTINE  PERFORMS  THE  RICCATI  TRANSFORM  FIRST 
C     SWEEP,  GENERATING  AND  STORING  ON  DISK  1  RMAT 
C     AND  SVEC  (FOR  EACH  MODE)  AT  EACH  FORWARD  STEP. 
C 

COMPLEX  RMAT(50,50),SVEC(50,100) 
COMPLEX  A(50,50),B(50,50),C(50,50),P(50,100) 
COMPLEX  D(50),SUM,DET 
REAL  COND,DMAG 

INTEGER  I,J,K,L,NA,NB,NC,PERND,BIND,IMX,INORM 
INTEGER  NBMX 

CHARACTER'-n  CHAR2 ,  CHARS,  CHAR4 
C 

C0MM0N/BLK2/PERND , BIND 
C0MM0N/BLK8/A,B,C,P 
C0MM0N/BLK9/CHAR2,  CHAR3 ,  CHAR4 
C 

C     LOADING  THEN  INVERTING  (B+A'-RMAT) 
C     USING  MINIMUM  MEMORY  SINGLE  MATRIX  TECHNIQUE 
C     DATE  1/29/80  FOR  THIS  CHANGE 
C 

C     SKIPPING  FIRST  A^'^R  (WHEN  A  =  0,  ZERO  R  MATRIX) 
IFd.EQ.  1)  THEN 

DO  10,  J  =  1,  NB 

DO  10,  K  =  1,  NB 

RMAT(J,K)  =  (0.  ,0.  ) 

10  CONTINUE 
C     RMAT  =  A-'--R 

ELSE 
C 

IF((CHAR2.EQ.  'l').OR.  (CHAR2.EQ.  'i'))  THEN 
WRITE(19,'^)  '  OLD  R  MATRIX' 
DO  11,  J  =  1,  5 

WRITE(19,1000)  (REAL(RMAT(J,K)),  K  =  1,  5) 

11  CONTINUE 
WRITE(19,*)  '  ' 
WRITE( 19,^0  '  A  MATRIX' 
DO  12,  J  =  1,  5 

WRITE(19,1000)  (REAL(A(J,K)),  K  =  1,  5) 

12  CONTINUE 
WRITE(19,''^)  '  ' 

ENDIF 


112 


DO  30,  K  =  1,  NB 

DO  20,  J  =  1,  NB 

D(J)  =  (0.  ,0.  ) 
DO  20,  L  =  1,  NA 

D(J)  =  D(J)  +  A(J,L)'''RMAT(L,K) 
20  CONTINUE 

DO  30,  J  =  1,  NB 

RMAT(J,K)  =  D(J) 
30         CONTINUE 
C 

IF((CHAR2.EQ.  'l').OR.  (CHAR2.EQ.  'i'))  THEN 
WRITE(19,''-)  '  NEW  R  MATRIX' 
DO  13,  J  =  1,  5 

WRITE(19,1000)  (REAL(RMAT(J,K)),  K  =  1,  5) 
13         CONTINUE 

WRITE(19,'V)  '  ' 
END  IF 


END  IF 


IF((CHAR2.EQ.  'l').OR.  (CHAR2.EQ.  'i'))  THEN 
WRITE(19,'"-)  '  B  MATRIX' 
DO  14,  J  =  1,  5 

WRITE(19,1000)  (REAL(B(J,K)),  K  =  1,  5) 

14  CONTINUE 
WRITE(19,''0  '  ' 

END  IF 
C 

C     RMAT  =  B  +  RMAT 
DO  40,  J  =  1,  NB 

DO  40,  K  =  1,  NB 

RMAT(J,K)  =  RMAT(J,K)  +  B(J,K) 
40    CONTINUE 
C 

IF((CHAR2.EQ.  'l').OR.  (CHAR2.EQ.  'i'))  THEN 
WRITE(19,--'-)  '  NEWEST  R  MATRIX' 
DO  15,  J  =  1,  NB 

WRITE(9,1000)  (REAL(RMAT(J,K)),  K  =  1,  NB) 
WRITE(19,1000)  (REAL(RMAT(J,K)),  K  =  1,  NB) 

15  CONTINUE 
WRITE(9,->)  '  ' 
WRITE(19,-'-)  '  • 

END  IF 
C     INVERTING  THE  MATRIX  (B  +  A-R) 

CALL  CSMINV( RMAT , NBMX , NB , DET , COND , INORM) 
C 

DMAG  =  CABS (DET) 
WRITE(^,")  I, NBMX, NB, DMAG, COND 
C 

C     COMPUTING  THE  NEW  S -VECTORS 
C 

C     SKIPPING  FIRST  A'-'-S  (A  =  0) 
IF  (I.EQ.  1)  THEN 

CONTINUE 
ELSE 

DO  70,  K  =  1,  PERND 
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DO  60,  J  =  1,  NB 

D(J)  =  (0.0,0.0) 
DO  60,  L  =  1,  NA 

D(J)  =  D(J)  +  A(J,L)'>SVEC(L,K) 
60  CONTINUE 

DO  70,  J  =  1,  NB 

P(J,K)  =  P(J,K)  -  D(J) 
70  CONTINUE 

END  IF 
C     FINAL  SVEC  MULTIPLICATION 

DO  100,  K  =  1,  PERND 
DO  90,  J  =  1,  NB 

D(J)  =  (0.0,0.0) 
DO  90,  L  =  1,  NB 

D(J)  =  D(J)  +  RMAT(J,L)*P(L,K) 
90  CONTINUE 

DO  100,  J  =  1,  NB 

SVEC(J,K)  =  D(J) 
100  CONTINUE 

C 

C     STORING  I+l  SVEC  ON  DISK  7 
DO  110  J  =  1,  NB 

WRITE(7)  (SVEC(J,K),  K  =  1,  PERND) 
110    CONTINUE 
C 

C     FINAL  RMAT  MULTIPLICATION 
IFd.EQ.  IMX)  RETURN 
DO  130,  J  =  1,  NB 

DO  120,  K  =  1,  NC 

D(K)  =  (0.0,0.0) 
DO  120,  L  =  1,  NB 

D(K)  =  D(K)  -  RMAT(J,L)"C(L,K) 
120        CONTINUE 

DO  130,  K  =  1,  NC 

RMAT(J,K)  =  D(K) 
130   CONTINUE 
C 

C     STORING  I+l  RMAT  ON  DISK  7 
DO  140,  J  =  1,  NB 
WRITE(7)  (RMAT(J,K),  K  =  1,  NC) 
140    CONTINUE 
C 

1000   FORMAT(10(E9.  3,1X)) 
C 

RETURN 
END 
C 
C 

SUBROUTINE  SWEEP( IMX, NABC,SURBC, LINE, CHARl ,U,BCOND,PSI ,ANS, CHARS) 
C     THIS  ROUTINE  PERFORMS  THE  RICCATI  TRANSFORM  BACKSWEEP 
C     FROM  I=IMX  TO  1=1,  RECALLING  RMAT  AND 
C     SVEC  FROM  DISK  7  AT  EACH  BACKSTEP  TO  FORM 
C     THE  NODE  VECTORS  PSI,  THEN  STORING  THESE  ON 
C     DISK  8  FOR  EACH  APPLIED  DIRICHLET  B.C. 

COMPLEX  RMAT(50, 100), PSI(50, 100), SVEC(50, 100), D( 100), TEMP 
COMPLEX  SURBC( 100) ,ANS( 100) ,LINE(50) ,ANSB(50) ,U( 100, 100) 
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COMPLEX  BCOND(IOO) 

REAL  ERRORP , ERRORD , TERRN , TERRD , ATSERR 
INTEGER  I,J,K,L,NABC(100,3),IMX,PERND,BIND 
CHARACTER'-n  CHARl,  CHARS 

C  OMMON/ B  LK2 / PERND , B I ND 


IF((CHAR5.EQ.  'M').0R.  (CHARS.  EQ.  'm'))  THEN 

RETURN 
END  IF 
C 

C     INITIAL  DISK  READ  AT  IMX  (R  =  0,  NOT  WRITTEN,  =>  S  IS  READ  FIRST) 
WRITE(*,'V)  '  ' 
DO  90,  I  =  IMX,  1,  -1 

WRITE ( ''Si OS 0)  (IMX-I+1),  IMX 
IF(I.EQ.  IMX)  THEN 

DO  10,  J  =  NABC(I,2),  1,  -1 
BACKSPACE  7 

READ(7)  (PSI(J,K),  K  =  1,  PERND) 
BACKSPACE  7 
10  CONTINUE 

DO  20,  J  =  1,  NABC(I,2) 
DO  IS,  K  =  1,  PERND 

U(IMX,K)  =  PSI(J,K) 
15  CONTINUE 

20  CONTINUE 

C  SUBSEQUENT  DISK  READS 

ELSE 
C  READ  R  MATRIX 

DO  30,  J  =  NABC(I,2),  1,  -1 
BACKSPACE  7 

READ(7)  (RMAT(J,K),  K  =  1,  NABC(I,3)) 
BACKSPACE  7 
30  CONTINUE 

C  READ  S  VECTOR 

DO  40,  J  =  NABC(I,2),  1,  -1 
BACKSPACE  7 

READ(7)  (SVEC(J,K),  K  =  1,  PERND) 
BACKSPACE  7 
40  CONTINUE 

C  MULTIPLY  RMAT  =  RMAT''-PSI 

DO  60,  J  =  1,  NABC(I,2) 
DO  SO,  K  =  1,  PERND 

D(K)  =  (0.0,0.0) 

DO  SO,  L  =  1,  NABC(I,3) 

D(K)  =  D(K)  +  RMAT(J,L)*PSI(L,K) 
SO  CONTINUE 

DO  60,  K  =  1,  PERND 

RMAT(J,K)  =  D(K) 
60  CONTINUE 

C  PSI  =  RMAT  +  SVEC 

DO  70,  J  =  1,  NABC(I,2) 
DO  70,  K  =  1,  PERND 

PSI(J,K)  =  RMAT(J,K)  +  SVEC(J,K) 
70  CONTINUE 

ANSB(I)  =  (0.0,0.0) 
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DO  80,  J  =  1,  NABC(I,2) 
DO  75,  K  =  1,  PERND 

IF((J.EQ.NABC(I,2)).0R.  (I.EQ.  1))  THEN 

U(I,K)  =  PSI(J,K) 
ELSEIF(J.  EQ.  1)  THEN 

U(2-"^IMX-I,K)  =  PSI(J,K) 
ELSEIF(J. EQ. (NABC(I,2)+l)/2)  THEN 

ANSB(I)  =  ANSB(I)  +  PSI( J,K)*BC0ND(K) 
ENDIF 
75  CONTINUE 

80  CONTINUE 

ENDIF 
90    CONTINUE 

WRITE(8,*)  '  ' 
C 

DO  98,  J  =  1,  PERND 

ANS(J)  =  (0.0,0.0) 
DO  94,  L  =  1,  PERND 

ANS(J)  =  ANS(J)  +  U(J,L)*BCOND(L) 
94         CONTINUE 
98    CONTINUE 
C 

TERRN  =  0.  0 
TERRD  =  0.  0 
DO  100,  I  =  1,  PERND 

ERRORP  =  (CABS(ANS(I)  -  SURBC( I) ) )**2 
ERRORD  =  (CABS(SURBC(I)))'>"-2 
TERRN  =  TERRN  +  ERRORP 
TERRD  =  TERRD  +  ERRORD 

WRITE(13,1020)  I,  BCOND(I),  SURBC(I),  ANS(I),  ERRORP 
100   CONTINUE 

WRITE(13,^'-)  '  ' 
ATSERR  =  (TERRN/TERRD)-'"'^0.5 
WRITE(13,1030)  ATSERR 

WRITE(''',-'0  'RMS  ERROR  (FOR  THE  PERIMETER)     =  ', ATSERR 
WRITE(13,-")  '  ' 
C 

IF((CHAR1.EQ.  'C').OR.  (CHARl.EQ.  'c'))  THEN 

RETURN 
ELSE 

TERRN  =  0.0 
TERRD  =  0.0 
DO  110,  I  =  2,  BIND+1 

ERRORP  =  (CABS(ANSB(I)  -  LINE(I)))**2 
ERRORD  =  (CABS(LINE(I)))*'>2 
TERRN  =  TERRN  +  ERRORP 
TERRD  =  TERRD  +  ERRORD 

WRITE(13,1025)  (I-l),  LINE(I),  ANSB(I),  ERRORP 
110        CONTINUE 

WRITE( 13,^0  '  ' 
ATSERR  =  (TERRN/TERRD)**o.5 
WRITE(13,1040)  ATSERR 

WRITE('^''0  'RMS  ERRROR  (FOR  BISECTION  SEGMENT)  =  ', ATSERR 
ENDIF 
C 

DO  120,  I  =  1,  PERND 
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WRITE(30,1060)  (U(I,J),  J  =  1,  PERND) 

120  CONTINUE 
C 

1020  F0RMAT(1X,I3,2X,3(E14.  8 , 1X,E14.  8,3X) ,F10.  6) 

1025  F0RiMAT(lX,I3,2X,4(E14.  8,2X),F10.  6) 

1030  F0RMAT(1X,'  RMS  ERROR  (FOR  THE  PERIMETER)     =  ' ,F12. 6) 

1040  F0RMAT(1X,'  RMS  ERROR  (FOR  BISECTION  SEGMENT)  =  ' ,F12. 6) 

1050  F0RMAT(1X,I3,'  TOTAL  BACKSWEEP  ROWS  OUT  OF  ',13,'  COMPLETED') 

1060  F0RMAT(1X,100(E8.  2,E8.  2,2X)) 


C 


RETURN 
END 


C 

CC 

C 

SUBROUTINE  CSMINV( A , NDIM , N , DETERM , COND , INORM) 
C 

C     INORM  -  FLAG  TO  NORMALIZE  COLUMNS  AND  ROWS  OF  MATRIX  A 
C     MATRIX  NORMALIZATION  BY  M.  A.  MORGAN 
C     APRIL  24,1978 
C 

C     A      -  MATRIX  TO  INVERT  -  INPUT/OUTPUT 

C     NDIM   -  -  INPUT 

C     N      -  -  INPUT 

C     DETERM  -  DETERMINATE  OF  A  -  OUTPUT 

C     COND   -  CONDITION  NUMBER  OF  A       -  OUTPUT 
C     INORM   -  INTEGER  NORMALIZATION  FLAG   -  INPUT 
C 
C 

COMPLEX  A(50, 50), PIV0T(50),AMAX,T, SWAP, DETERM, U 
INTEGER  I,J,K,L,IPIVOT(50),INDEX(50,2),IROW,ICOLUM,L1,JROW 
INTEGER  JCOLUM,N, INORM 

REAL  TEMP,ALPHA(50),COL(50),ROW(50),AJK,SUMAXA,SUMROW,SUMAXI 
C 

IF(NDIM. GT.  50)  THEN 

WRITE(",'--)  '  ERROR  IN  INVERTION  CALL...  DIMENSION  >  50  ' 
STOP 
END  IF 
IF(N.  GT.  NDIM)  THEN 

WRITE(^'S^O  '  ERROR  IN  INVERTION  CALL...  N  >  MAX  DIM.  ' 
STOP 
END  IF 

IF( INORM. NE. 1)  GO  TO  7 
DO  3  K  =  1,  N 

COL(K)  =0.0 
DO  1  J  =  1,N 

AJK  =  CABS(A(J,K)) 

IF(AJK. GT. COL(K))  COL(K)  =  AJK 

1  CONTINUE 

DO  2  J  =  1,  N 

2  A(J,K)  =  A(J,K)/COL(K) 

3  CONTINUE 

C     ROW  NORMALIZING 
DO  6  J  =  1,  N 

ROW(J)  =0.0 
DO  4  K  =  1,  N 
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AJK  =  CABS(A(J,K)) 

IF(AJK. GT. ROW(J))  ROW(J)  =  AJK 

4  CONTINUE 

DO  5  K  =  1,  N 

5  A(J,K)  =  A(J,K)/ROW(J) 

6  CONTINUE 

7  CONTINUE 

DETERM  =  CMPLX(  1.0,0.0) 
SUMAXA  =  0.  0 
DO  20  J  =  1,  N 

ALPHA(J)  =0.0 
SUMROW  =  0.  0 
DO  10  I  =  1,  N 

ALPHA(J)  =  ALPHA(J)  +  A(J,I)*  CONJG( A( J, I) ) 
10  SUMROW  =  SUMROW  +  CABS  (A(  J, I)) 

ALPHA(J)  =  SQRT(ALPHA(J)) 
IF(SUMROW.GT. SUMAXA)  SUMAXA  =  SUMROW 
20    IPIVOT(J)  =  0 

DO  600  I  =  1,  N 
C 

c 

AMAX  =  CMPLX(0.  0,0.  0) 
DO  105  J  =  1,  N 

IF  (IPIVOT(J)-l)  60,  105,  60 
60  DO  100  K  =  1,  N 

IF  (IPIVOT(K)-l)  80,  100,  740 
80  TEMP  =  AMAX'^CONJG(AMAX)  -  A(  J,K)'>CONJG(  A(  J,K)  ) 

IF(TEMP)85,85,100 
85  IROW  =  J 

ICOLUM  =  K 
AMAX  =  A(J,K) 
100  CONTINUE 

105        CONTINUE 

IPIVOT( ICOLUM)  =  IPIVOTC ICOLUM)  +  1 
C 

c 

IF  (IROW-ICOLUM)  140,  260,  140 
140  DETERM  =  -DETERM 

DO  200  L  =  1,  N 

SWAP  =  A(IROW,L) 
A(IROW,L)  =  A( ICOLUM, L) 
200  A( ICOLUM, L)  =  SWAP 

SWAP  =  ALPHA(IROW) 

ALPHA(IROW)  =  ALPHA( ICOLUM) 

ALPHA( ICOLUM)  =  SWAP 
260  INDEX(I,1)  =  IROW 

INDEX(I,2)  =  ICOLUM 

PIVOT(I)  =  A( ICOLUM, ICOLUM) 

U  =  PIVOT(I) 

DETERM  =  DETERM-U 

DETERM  =  DETERM/ALPHA( ICOLUM) 

TEMP  =  PIVOT(I)''-CONJG(PIVOT(I)) 

IF(TEMP)330, 720,330 
C 

c 

330  AdCOLUM,  ICOLUM)   =  CMPLX(  1.  0  ,0.  0) 
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DO  350  L  =  1,  N 

U  =  PIVOT(I) 
350  A(ICOLUM,L)  =  A( ICOLUM,L)/U 

C 

c 

380  DO  550  LI  =  1,  N 

IF(Ll-ICOLUM)  400,  550,  400 
400  T  =  A(L1,IC0LUM) 

A(L1,IC0LUM)  =  CMPLX(0.  0,0.  0) 
DO  450  L  =  1,  N 

U  =  A(ICOLUM,L) 
450  A(L1,L)  =  A(L1,L)  -  U*T 

550        CONTINUE 
600   CONTINUE 
C 
C 
620   DO  710  I  =  1,  N 

L  =  N  +  1  -  I 

IF  (INDEX(L,1)  -  INDEX(L,2))  630,  710,  630 
630  JROW  =  INDEX(L,1) 

JCOLUM  =  INDEX(L,2) 
DO  705  K  =  1,  N 

SWAP  =  A(K,JROW) 
A(K,JROW)  =  A(K, JCOLUM) 
A(K, JCOLUM)  =  SWAP 
705  CONTINUE 

710   CONTINUE 

SUMAXI  =  0. 0 
DO  910  I  =  1,  N 

SUMROW  =  0. 0 
DO  900  J  =  1,  N 
900        SUMROW  =  SUMROW  +  CABS(A(I,J)) 

IF(SUMROW.GT. SUMAXI)  SUMAXI  =  SUMROW 
910   CONTINUE 

COND  =  SUMAXA-SUMAXI 
IF(INORM.NE. 1)  GO  TO  955 
DO  950  K  =  1,  N 

DO  950  J  =  1,  N 
950        A(J,K)  =  A(J,K)/(ROW(K)^-COL(J)) 
955   CONTINUE 

RETURN 
720   WRITE(''S730) 

730   FORMAT(  '  MATRIX  IS  SINGULAR') 
740   RETURN 
END 

SUBROUTINE  SAVE( BCOND , ANS , U , OFFSET , PERND , CHAR5 , DPER , K , XORG , YORG , 
CNRES,MRES,LBIAS,GBIAS,MAXD) 
C 

C     THIS  SUBROUTINE  SAVES  THE  ESSENCE  OF  THE  FINITE  ELEMENT  PROBLEM 
C     TO  A  DATA  FILE  CALLED  "  F3. DAT  ".   tHIS  DATA  IS  NECESSARY  TO 
C     SOLVE  THE  FIELD  FEEDBACK  FORMULATION,  (F3). 
C 

COMPLEX  BCOND( 100) , ANS( 100) ,U( 100 , 100) 

REAL  OFFSET, MESH(0: 200,5) , DPER, K, XORG, YORG, MRES ,MAXD 

INTEGER  PERND , I , J , NRES , LBI AS , GBI AS 
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CHARACTER*!  CHARS 

COMMON/BLKl/MESH 

IF(  (CHARS.  EQ.  'M').0R.  (CHAR5.EQ.  'm'))  THEN 

RETURN 
ENDIF 

WRITE(40,'^)  PERND 

WRITE ( 40, '>)  OFFSET 

WRITE ( 40, ''0  DPER 

WRITE(40,''O  K 

WRITE ( 40, -'0  XORG 

WRITE (40,''')  YORG 

WRITE(40,''-)  NRES 

WRITE(40,''-)  MRES 

WRITE(40,''O  LBIAS 

WRITE ( 40, ''0  GBIAS 

WRITE (40,''-)  MAXD 

DO  10,  I  =  1,  4 

DO  10,  J  =  1,  PERND 
WRITE(40,''O  MESH(J,I) 


10 

CONTINUE 

c 

DO  20,  J  =  1,  PERND 

WRITE(40,-'O  BCOND(J) 

20 

CONTINUE 

C 

DO  30,  J  =  1,  PERND 

WRITE(40,-'O  ANS(J) 

30 

CONTINUE 

C 

DO  40,  1=1,  PERND 

DO  40,  J  =  1,  PERND 

WRITE(40,''O  U(I,J) 

40 

CONTINUE     ---^ 

C 

RETURN 

END 

c 

c 
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APPENDIX  D.     VARINT  CONVERGENCE  PROGRAM 

C 

C       TEST  OF  VARINT  CONVERGENCE 
C 

COMPLEX  F(3,3),  ALPHA,  BETA,  EXACT,  CENTER,  LEFT,  RIGHT,  TOP 
COMPLEX  BOTTOM,  SUM,  CALC 
REAL  X(3),  Y(3),  D,  AREA,  KR,  PI,  ERROR 
INTEGER  I 

OPEN  (UNIT  =  1,  FILE  =  'C:  MSFORT  TEST.  DAT',  STATUS  =  'UNKNOWN') 
ALPHA  =  (1.0,0.0) 
BETA  =  (1.0,0.0) 
PI  =  4.  O^^ATANd.  0) 
C 

DO  5,  1=1,  100 

D  =  2. 0"PI"FLOAT(I)/100.0 
KR  =  REAL(CSQRT(BETA)) 
X(l)  =  0.  0 
Y(l)  =  0.  0 
X(2)  =  D 
Y(2)  =0.0 
X(3)  =  0.  0 
Y(3)  =  D 
C 

CALL  VARINT(X,Y,F, ALPHA, BETA, AREA) 
C 

EXACT  =  CMPLX(COS(KR"-0.0),SIN(KR-'-0.0)) 
CENTER  =  4.  0---F(l,l) 

RIGHT  =  -2.  0^'-F(l,2)^'^CMPLX(COS(KR-"-Y(2)),SIN(KR"-Y(2))) 
TOP  =  -2.  0-->F(1,2)^>CMPLX(COS(KR'-Y(3)),SIN(KR-->Y(3))) 
LEFT  =  -2.  0-''-F(l,2)^'^CMPLX(C0S(KR^'^Y(2)),SIN(KR^--Y(2))) 
BOTTOM  =  -2.  0*F(l,2)'^CMPLX(COS(-KR'-Y(3)),SIN(-KR-'-Y(3))) 
SUM  =  TOP  +  BOTTOM  +  LEFT  +  RIGHT 
CALC  =  SUM/CENTER 
ERROR  =  (CALC -EXACT) /EXACT 
WRITE( 1,1000)  I, D, EXACT, CALC, ERROR 
5     CONTINUE 
C 

CLOSE(l) 
C 

1000  F0RMAT(1X,I3,1X,F8.5,1X,2(F8.5,1X,F8.5,1X),E13.6) 
C 

STOP 
END 
C 
C 

SUBROUTINE  VARINT( X , Y , F , ALPHA , BETA , AREA ) 
C 

C     GENERATING  VARIATIONAL  FINITE  ELEMENT  AREA  INTEGRATIONS  OF  THE 
C     LINEAR  BASIS  FUNCTION  LAGRANGIAN  FOR  THE  HELMHOLTZ  EQUATION. 
C     THESE  ARE  RETURNED  IN  F(3,3).  X(3)  AND  Y(3)  ARE  THE  WAVENUMBER 
C     NORMALIZED  CARTESIAN  COORDINATES  OF  THE  TRIANGLE  VERTICES. 
C     X  =  Ko-'^x,  Y  =  Ko-''--y    ALPHA  AND  BETA  ARE  COMPLEX 


121 


10 

c 


MATERIAL  PARAMETERS  WITHIN  THE  ELEMENT. 

FOR  TM  INCIDENCE:  ALPHA  =  1/ur;  BETA  =  er 
FOR  TE  INCIDENCE:  ALPHA  =  1/er;  BETA  =  ur 

OUTPUTS  :    F(3,3)  -  FINITE  ELEMENT  AREA  INTEGRATION 
AREA   -  AREA  OF  A  TRIANGLE 


COMPLEX  ALPHA,  BETA,  B12,  F(3,3) 
REAL  X(3),  Y(3),  T(3,3),  AREA,  DET 
INTEGER  L,  K 

DET  =  ABS(X(2)*Y(3)  +  X(3)''-Y(l)  +  X(1)*Y(2) 
CX(1)''-Y(3)  -  X(2)-''^Y(1)) 
AREA  =  ABS(0.5''^DET) 
B12  =  BETA/12. 

Y(3))/DET 
Y(1))/DET 
Y(2))/DET 
X(2))/DET 
X(3))/DET 
X(1))/DET 

X(3)''^Y(2))/DET 
X(1)'>Y(3))/DET 


X(3)*Y(2)  - 


T(l,l) 
T(l,2) 
T(l,3) 
T(2,l) 
T(2,2) 
T(2,3) 
T(3,l) 
T(3,2) 
T(3,3) 
DO  10, 
DO 


=  (Y(2) 
=  (Y(3) 
=  (Yd) 
=  (X(3) 
=  (X(l) 
=  (X(2) 
=  (X(2)'''Y(3) 
=  (X(3)-'VY(1) 
=  (X(l)''-Y(2) 
K  =  1,  3 
10  L  =  1,  3 

F(K,L)  =  ALPHA->(T(1,K)'>T(1,L)  +  T(2,K)*T(2  ,L) ) 
IF(K. EQ. L)  F(K,L)  =  F(K,L)  -  B12 
F(K,L)  =  AREA--'^F(K,L) 

RETURN 
END 


X(2)''-Y(1))/DET 


B12 
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APPENDIX  E.     FIELD  FEEDBACK  PROGRAM 

c 

C     FIELD  FEEDBACK  FORMULATION  PROGRAM 

C     WRITTEN  BY  T.  B.  WELCH 

C 

C     w/  PROGRAMMING  IDEAS  FROM  PROF  M.  A.  MORGAN 

C 

C     BCOND     -  BOUNDARY  CONDITIONS 

C     OFFSET    -  OFFSET  IN  WAVELENGTHS  (PERIMETER  TO  BOUNDARY) 

C     ANS       -  CALCULATED  PS I  VALUES  ON  PERIMETER 

C     DPER      -  DESIRED  PERCENT  ERROR  SCALE  FACTOR  FOR  GREEN'S 

C  FUNCTION  INTEGRAL  PATCH  STEPPING 

C     LBIAS     -  BIAS  THAT  IS  ADDED  TO  THE  GFI  STEP  FOR  <  1. 0 

C     GBIAS      -  BIAS  THAT  IS  ADDED  TO  THE  GFI  STEP  FOR  >  1.  0 

C     MAXD      -  MAXIMUM  DISTANCE  BEYOND  WHICH  NO  CONTRIBUTION  IS  MADE 

C  TO  THE  GFI 

C     K         -  WAVENUMBER 

C     XORG      -  X  ORIGIN 

C     YORG      -  Y  ORIGIN 

C     NRES      -  NUMBER  OF  EVENLY  SPACED  POINTS  DESIRED  FOR  THE  FAR 

C  FIELD  CALCULATIONS  (360/NRES  =  ANGULAR  RESOLUTION) 

C     U         -  MATRIX  THAT  RELATES  EACH  BOUNDARY  NODE  VALUE  TO  THE 

C  UNKNOWN  PERIMETER  NODE  VALUE.   MULTIPLY  U  BY  A  DRIVING 

C  VECTOR  (ON  BOUNDARY)  TO  FIND  PERIMETER  VALUES. 

C     PERND     -  NUMBER  OF  PERIMETER  NODES 

C     MESH      -  GEOMETRY  ARRAY  CONTAINING: 

C  1  -  X  POSITION  OF  PERIMETER  NODES 

C  2  -  Y  POSITION  OF  PERIMETER  NODES 

C  3  -  X  UNIT  NORMAL  OF  PERIMETER  NODES 

C  4  -  Y  UNIT  NORMAL  OF  PERIMETER  NODES 

C  5  -  X  POSITION  OF  GEOMETRIC  CONTOUR  NODES 

C  6  -  Y  POSITION  OF  GEOMETRIC  CONTOUR  NODES 

C  7  -  X  POSITION  OF  BOUNDARY  NODES 

C  8  -  Y  POSITION  OF  BOUNDARY  NODES 

C     T         -  MATRIX  THAT  RELATES  PERIMETER  VALUES  BACK  OUT 

C  TO  THE  BOUNDARY  VIA  A  GREEN'S  FUNCTION  INTEGRAL 

C     CNVEC     -  VECTOR  OF  SCATTERED  FIELD  BACK  ONTO  THE  BOUNDARY 

C     MRES      -  MESH  RESOLUTION 

C 

C 

COMPLEX  BCOND( 100) ,ANS( 100) ,U( 100 , 100) ,T( 100, 100) ,CNVEC( 100) 

REAL  OFFSET, MESH( 100, 8), DPER, K, XORG, YORG, MRES, MAXD 

INTEGER  PERND, I, J, NRES, LB I AS, GBIAS 
C 

OPEN  (UNIT  =  40,  FILE  =  'C:  MSFORT  F3.  DAT' ,STATUS=' UNKNOWN ' ) 

OPEN  (UNIT  =  50,  FILE  =  'C:  MSFORT  FFPAT.  DAT' ,STATUS=' UNKNOWN ' ) 
C 

CALL  INPUT( BCOND, ANS, U, OFFSET, PERND, MESH, DPER, K, NRES, XORG, YORG, 
CMRES, LBIAS, GBIAS, MAXD) 

CALL  TMAT(U, PERND, MESH, T, OFFSET, DPER, BCOND, MRES, LBIAS, GBIAS, MAXD) 

CALL  CNSOLV(T, BCOND, CNVEC, PERND) 
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CALL  FFLD(CNVEC,PERND, MESH, U, OFFSET, K,NRES,XORG,YORG) 

CLOSE(40) 
CL0SE(50) 

STOP 
END 


SUBROUTINE  INPUT( BCOND , ANS , U , OFFSET , PERND , MESH , DPER , K , NRES , XORG , 
CYORG , MRES , LBI AS , GBIAS , MAXD) 

THIS  SUBROUTINE  READS  THE  FINITE  ELEMENT  PROBLEM  DATA  FROM 
THE  DATA  FILE  CALLED  "  F3. DAT  ".   THIS  DATA  IS  NECESSARY  TO 
SOLVE  THE  FIELD  FEEDBACK  FORMULATION,  (F3). 

COMPLEX  BCOND( 100) ,ANS( 100) ,U( 100, 100) 

REAL  OFFSET, MESH( 100, 8), DPER, K, XORG, YORG, MRES, MAXD 

INTEGER  PERND, I, J, NRES, LBIAS, GBIAS 

WRITE(*,-v)  '  READING  INPUT  DATA  ' 


READ( 
READ( 
READ( 
READ( 
READ( 
READ( 
READ( 
READ( 
READ( 
READ( 
READ( 


40 
40 
40 
40 
40 
40 
40 
40 
40 
40 
40 


WRITE(6 
WRITE (6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE (6 
WRITE(6 
WRITE (6 
WRITE(6 
WRITE(6 


*)   PERND 
•'')    OFFSET 
■''')    DPER 

*)  XORG 

''-)  YORG 

-'0  NRES 

-'-)  MRES 

''0  LBIAS 

^0  GBIAS 

*)  MAXD 

*)  '  NUMBER  OF  PERIMETER  NODES 

*)  'BOUNDARY  CONTOUR  OFFSET 

*)  '  DESIRED  GFI  SCALE  FACTOR 

-'■)  '  GFI  STEP  BIAS  FOR  <  1.  0 

''0  '  GFI  STEP  BIAS  FOR  >  1.0 

'■)  '  MAX  DIST  >  NO  CONTRIBUTION  TO  GFI 

''0  '  WAVENUMBER 

'■')  '  X  ORIGIN 

*)  '  Y  ORIGIN 

''0  '  NUMBER  OF  NODES  FOR  SIGMA 

*)  '  REQUESTED  MESH  RESOLUTION 


, PERND 

, OFFSET 

,DPER 

, LBIAS 

, GBIAS 

,MAXD 

,K 

,XORG 

,YORG 

,NRES 

,MRES 


DO  10,  I  =  1,  4 

DO  10,  J  =  1,  PERND 

READ(40,*)  MESH(J,I) 


10 

c 

20 
C 

CONTINUE 

DO  20,  J  =  1,  PERND 

READ(40,*)  BCOND(J) 
CONTINUE 

DO  30,  J  =  1,  PERND 
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30 
C 


40 
C 


50 
C 


MESH(  I,  3)  •'^OFFSET/2.  0 
MESH(I,4)^^0FFSET/2.0 
MESH(I,3)''^0FFSET 
ME  SH(  1, 4) ''^OFFSET 


READ(40,''O  ANS(J) 

CONTINUE 

DO  40,  1=1,  PERND 

DO  40,  J  =  1,  PERND 

READ(40,^'O  U(I,J) 
CONTINUE 

DO  50,  1=1,  PERND 

MESH(I,5)  =  MESH(I,1)  + 

MESH(I,6)  =  MESH(I,2)  + 

MESH(I,7)  =  MESH(I,1)  + 

MESH(I,8)  =  MESH(I,2)  + 

CONTINUE 

RETURN 

END 


SUBROUTINE  TMAT( U , PERND , MESH , T , OFFSET , DPER , BCOND , MRES , LEI AS , GB IAS , 
CMAXD) 

THIS  SUBROUTINE  CALCULATES  THE  GREEN'S  FUNCTION  INTEGRAL 

(FOR  A  SINGLE  BASIS  FUNCTION  BOUNDARY  CONDITION)  GEOMETRIC 

PERIMETER  WITH  RESPECT  TO  EACH  OF  THE  OFFSET  BOUNDARY  NODES. 

THE  INTEGRATION  IS  REPEATED  UNTIL  EACH  BOUNDARY  CONDITION  HAS  BEEN 

INDIVIDUALLY  APPLIED  AND  INTEGRATED  WITH  RESPECT  TO  EACH  OFFSET 

BOUNDARY  NODE.   THESE  VALUES  ARE  RETURNED  IN  THE  "  T  "  MATRIX 


FOR  USE  IN  THE  FIELD  FEEDBACK  FORMULATION. 
ORGANIZED,  T  m,n. 


THE  MATRIX  IS 


COMPLEX  U( 100, 100) ,T( 100 , 100) ,PVEC( 100) ,PDVEC( 100) 

COMPLEX  J,HORP1,HORP2,H1RP1,SUM,TEMP(100),BCOND(100) 

COMPLEX  H1RP2,PSI,PSIRP, INTEGRAL, DPSIC 

REAL  RMRP,N0RM11,N0RM12, DOT, DPER, DISTM,MAXD 

REAL  OFFSET, MESHdOO, 8), DIST,R,DZ,DL, MRES 

INTEGER  I , M , N , NN , STEP , FN , SN , PERND , MM , STEPMX , STEPMN , LB IAS , GBIAS 

OPEN  (UNIT  =   2,  FILE  =  'C:  MSFORT  TMAT. DAT' , STATUS  =  'UNKNOWN') 

J  =  (0.0,1.0) 

STEPMX  =  INT(DPER^'-(-48.0^^0FFSET  +  17.2)  +  LBIAS) 

STEPMN  =  INT(DPER  +  GBIAS) 


WRITE(''-,''0  '  LOADING  T  MATRIX  ' 

WRITE(''-,*)  •  MAXIMUM  STEP  =  ', STEPMX,',  MINIMUM  STEP  =  ', STEPMN 
WRITE(",^'-)  •  MAXIMUM  DISTANCE  FOR  ANY  CONTRIBUTION  =  '  ,MAXD 
DO  40,  M  =  1,  PERND 

DO  5,  1=1,  PERND 

IF(M.  EQ.  I)  THEN 

PVEC(I)  =  (1.0  +  U(I,M))/2.0 
PDVEC(I)  =  (1.0  -  U(I,M))/OFFSET 
ELSE 

PVEC(I)  =  U(I,M)/2.0 
PDVEC(I)  =   -U( I, M) /OFFSET 
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END  IF 
5  CONTINUE 

C 

DO  30,  N  =  1,  PERND 

WRITE('^IOOO)  M,  N,  PERND 
SUM  =  (0.0,0.0) 
DO  20,  NN  =  1,  PERND 
FN  =  NN 
IF(NN.EQ.  PERND)  THEN 

SN  =  1 
ELSE 

SN  =  NN  +  1 
ENDIF 

DIST  =  SQRT((MESH(FN,5)-MESH(SN,5))**2+(MESH(FN,6) 
C  -MESH(SN,6))'^'>2) 

INTEGRAL  =  (0.0,0.0) 
C 

DISTM  =  SQRT((MESH(N,7)-MESH(FN,5))'^'>2+(MESH(N,8)-MESH(FN,6))**2) 
C 

R  =  MESH(SN,5)  -  MESH(FN,5) 
DZ  =  MESH(SN,6)  -  MESH(FN,6) 
DL  =  SQRT(R-->--'-2  +  DZ-''"'>-2) 
NORMll  =  -DZ/DL 
NORM 12  =  R/DL 
C 

IF(DISTM. GT.  MAXD)  THEN 

GOTO  20 
ELSEIF(DISTM.  LE.  1.  0)  THEN 

STEP  =  STEPMX 
ELSE 

STEP  =  STEPMN 
ENDIF 

IF(STEP.  LT.  1)  STEP  =  1 
C 

DO  10,  I  =  1,  STEP+1 
IFd.EQ.  1)  THEN 

RMRP  =  SQRT((MESH(N,7)-(MESH(FN,5)+0.25'>(MESH(SN,5)- 
C  MESH(FN,5))/FL0AT(STEP)))""-2+(MESH(N,8)-(MESH(FN,6)  + 

C  0.  25->(MESH(SN,6)-MESH(FN,6))/FL0AT(STEP)))"-"2) 

DPSIC  =  PDVEC(FN)  +  0. 25"( PDVEC( SN) -PDVEC(FN) )/STEP 
PSIRP  =  PVEC(FN)  +  0.  25''^(PVEC(SN)-PVEC(FN))/STEP 
DOT  =  (N0RM11^--(MESH(N,7)  -  (MESH(FN,5)+0.  25"(MESH(  SN,5) - 
C  MESH(FN,5))/STEP))+(N0RM12''-(MESH(N,8)-(MESH(FN,6)+0.  25-^ 

C  (MESH(SN,6)-MESH(FN,6))/STEP))))/RMRP 

ELSEIF( I.  EQ.  STEP+1)  THEN 

RMRP  =  SQRT((MESH(N,7)-(MESH(FN,5)+(FLOAT(STEP)-0.25)"( 
C  MESH(  SN ,  5  )  -MESH(  FN ,  5  )  )  /FLOAT(  STEP)  )  )-'"'-2+(  MESH(  N ,  8 )  - ( MESH 

C  .        (FN,6)+(FLOAT(STEP)-0.25)''^(MESH(SN,6)-MESH(FN,6))/FLOAT 
C  (STEP))) -'"'-2) 

DPSIC=PDVEC(FN)+(FL0AT(STEP)-0.25)*(PDVEC(SN)-PDVEC(FN)) 
C  /(STEP) 

PSIRP  =  PVEC(FN)+(FLOAT(STEP)-0.25)*(PVEC(SN)-PVEC(FN))/ 
C  (STEP) 

DOT  =  (NORM11^'^(MESH(N,7)-(MESH(FN,5)+(FLOAT(STEP)-0.  25)" 
C  (MESH(SN,5)-MESH(FN,5))/STEP))+(N0RM12'>(MESH(N,8)-(MESH 

C  (FN,6)+(FL0AT(STEP)-0.  25)''-(MESH(  SN,6) -MESH(FN,6)  )/STEP) 
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C  )))/RMRP 

ELSE 

RMRP  =  SQRT((MESH(N,7)-(MESH(FN,5)+FLOAT(I-l)^nMESH(SN, 
C  5)-MESH(FN,5))/FL0AT(STEP)))""2+(MESH(N,8)-(MESH(FN,6)+ 

C  FLOAT(  I  - 1  )'>(  MESH(  SN ,  6 )  -MESH(  FN ,  6 )  )  /FLOAT(  STEP)  )  )'>--'^2 ) 

DPSIC=PDVEC(FN)+(I-1)^'-(PDVEC(SN)-PDVEC(FN))/(STEP) 
PSIRP  =  PVEC(FN)+(I-1)''^(PVEC(SN)-PVEC(FN))/(STEP) 
DOT  =  (N0RM11"(MESH(N,7)-(MESH(FN,5)+(I-1)^KMESH(SN,5)- 
C  MESH(FN,5))/STEP))+(N0RM12-''-(MESH(N,8)-(MESH(FN,6)+(I-l)* 

C  ( MESH( SN , 6 ) -MESH( FN , 6 ) ) /STEP) ) ) ) /RMRP 

END  IF 

CALL  HAN 1(CMPLX( RMRP, 0. 0) ,H0RP1 ,H1RP1) 
IF(  (  I .  EQ.  1 )  .  OR.  ( I .  EQ.  STEP+1 )  )  THEN 

PSI  =  (J/4.  0)-'KHORPl^>-DPSIC  -  PSIRP'''D0T''^HlRPl)/2 
ELSE 

PSI  =  (J/4.  0)^--(HORPl''^DPSIC  -  PSIRP^'-DOT^'^HIRPI) 
END  IF 

INTEGRAL  =  INTEGRAL  +  PSI*DIST/STEP 
C 
10  CONTINUE 

SUM  =  SUM  +  INTEGRAL 
20  CONTINUE 

T(N,M)  =  SUM 
30         CONTINUE 
40    CONTINUE 
C 

DO  55,  1=1,  PERND 

TEMP(I)  =  (0. 0,0. 0) 
DO  50,  M  =  1,  PERND 

TEMP(I)  =  TEMP(I)  +  T(I,M)^--BCOND(M) 
50         CONTINUE 

WRITE(2,--)  (T(I,MM),  MM  =  1 ,  PERND) 
55    CONTINUE 
C 

DO  60,  1=1,  PERND 

WRITE (2,^0  I,  TEMP(I) 
60    CONTINUE 
C 

1000  FORMAT( IX ,' COLUMN  ',13,',  ROW  ',13,'  OUT  OF  ',13) 
C 

CL0SE(2) 
C 

RETURN 
END 
C 
C 

SUBROUTINE  CNSOLV( T , BCOND , CNVEC , PERND ) 
C 

C     THIS  SUBROUTINE  CALCULATES  THE  C  VECTOR  BY  SOLVING: 
C 

C  -1 

C     Cn  =  [I  -  T]   *  BOUNDARY  CONDITIONS  (INCIDENT  FIELDS) 
C 

COMPLEX  BCOND(100),T(100,100),TEMP(100),CNVEC(100),  DETERM 
REAL  COND,  DMAG 
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INTEGER  PERND,I,J,K,L,INORM,NMAX 

I NORM  =  0 
NMAX  =100 


DO  10,  I  =  1,  PERND 

DO  10,  J  =  1,  PERND 
IF(I.EQ.  J)  THEN 

T(I,J)  =  1  -  T(I,J) 
ELSE 

T(I,J)  =  -T(I,J) 
END  IF 
10    CONTINUE 

WRITE('^*)  'INVERTING  THE  [  I  -  T]  MATRIX 
C 

CALL  CSMINV( T , NMAX , PERND , DETERM , COND , INORM) 
C 

DMAG  =  CABS (DETERM) 

WRITE('>,-'-)  '  DETERMINANT      =  ',  DMAG 
WRITE('^''0  '  CONDITION  NUMBER  =  ',  COND 

WRITE('>,''0  '  MULTIPLING  MATRICES  TO  FORM  THE  SCATTERED  FIELDS  ' 
C 

DO  30,  1=1,  PERND 

CNVEC(I)  =  (0.  0,0.  0) 
DO  30,  J  =  1,  PERND 

CNVEC(I)  =  CNVEC(I)  +  T( I , J)*BCOND( J) 
30    CONTINUE 
C 

RETURN 
END 
C 
C 

SUBROUTINE  FFLD(CNVEC, PERND, MESH, U, OFFSET, K,NRES,XORG,YORG) 
C 

C     THIS  SUBROUTINE  CALCULATES  THE  FAR  FIELDS  DUE  TO  THE  OFFSET 
C     BOUNDARY  SCATTERED  FIELDS  AND  THE  PERIMETER  SCATTERED  FIELDS. 
C     ADDITIONAL  GREEN'S  FUNCTION  INTEGRALS  ARE  ACCOMPLISHED. 
C 

COMPLEX  CNVEC( 100) ,U( 100 , 100) , J,PSISP( 100) ,TEMP,PSI ,DPSI 
COMPLEX  DPSISP( 100) , INTEGRAL, DPSIM( 100 , 100) ,PSIM( 100, 100) 
REAL  MESH( 100, 8), OFFSET, K, DOT, DOTl, PI, ARES, XORG,YORG,DIST, SIGMA 
REAL  R,DZ,DL,N0RM11,N0RM12 
INTEGER  PERND,I,M,N,L,NRES,FN,SN,STEP,II 
C 

WRITE('^''0  '  CALCULATING  THE  SCATTERED  FIELDS' 
J  =  (0.0,1.0) 
PI  =  4.  0->ATAN(l.  0) 
ARES  =  2.  0^--PI/FLOAT(NRES) 
STEP  =  5 
C 

DO  5,  M  =  1,  PERND 

DO  5,  1=1,  PERND 

IF(M.  EQ.  I)  THEN 

PSIM(I,M)  =  (1.0  +  U(I,M))/2.0 
DPSIM(I,M)  =  (1.0  -  U(I,M))/OFFSET 
ELSE 
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PSIM(I,M)  =  U(I,M)/2.0 
DPSIMCI,M)  =   -U(I,M)/OFFSET 
ENDIF 


5 
C 

CONTINUE 

DO  10,  1=1,  PERND 

PSISP(I)  =  (0.0,0.  0) 

DPSISP(I)  =  (0.0,0.0) 

DO  10,  L  =  1,  PERND 

PSISP(I)  =  PSISP(I)  + 

DPSISP(I)  =  DPSISP(I) 

10 

CONTINUE 

C 

PSIM(I,L)'>CNVEC(L) 
+  DPSIM(I,L)'>CNVEC(L) 


DO  30,  1=0,  NRES-1 

WRITE(6,1000)  I+l,  NRES 
INTEGRAL  =  (0.0,0.0) 
DO  20,  N  =  1,  PERND 
FN  =  N 
IF(N.EQ.  PERND)  THEN 

SN  =  1 
ELSE 

SN  =  N  +  1 
ENDIF 

R  =  MESH(SN,5)  -  MESH(FN,5) 
DZ  =  MESH(SN,6)  -  MESH(FN,6) 
DL  =  SQRT(R-''"'^2  +  DZ''">2) 
NORMll  =  -DZ/DL 
NORM 12  =  R/DL 
DO  15,  II  =  1,  STEP+1 

DIST  =  SQRT((MESH(FN,5)-MESH(SN,5))*''-2+(MESH(FN,6)- 
C  MESH(SN,6))-''^^'^2) 

IFdI.EQ.  1)  THEN 

PSI  =  PSISP(FN)+0.25*(PSISP(SN)-PSISP(FN))/ 
C  FLOAT(STEP) 

DPSI  =  DPSISP(FN)+0.  25->(DPSISP(SN)-DPSISP 
C  (FN))/FLOAT(STEP) 

DOT  =  NORMll-"-SIN(I'''ARES)  +  NORM  1 2'''C0S ( I^'^ ARES) 
DOTl  =  (MESH(FN,5)+0.  25^'^(MESH(SN,5)-MESH(FN, 
C  5))/FL0AT(STEP)-X0RG)-'^SIN(I''^ARES)  +  (MESH(FN,6)  + 

C  0.  25-"(MESH(SN,6)-MESH(FN,6))/FL0AT(STEP)- 

C  YORG)-'^COS(I^^ARES) 

ELSEIF( II.  EQ.  STEP+1)  THEN 

PSI  =  PSISP(FN)+(FLOAT(STEP)-0.  25)'^(PSISP(SN)- 
C  PSISP(FN))/FLOAT(STEP) 

DPSI  =  DPSISP(FN)+(FLOAT(STEP)-0.  25)*(DPSISP(SN)- 
C  DPSISP(FN))/FLOAT(STEP) 

DOT  =  NORM  11--'^SIN(I''^ ARES)  +  N0RM12*C0S(  I^^areS) 
DOTl  =  (MESH(FN,5)+(FL0AT(STEP)-0. 25)"(MESH(SN,5)- 
C  MESH(FN,5))/FL0AT(STEP)-X0RG)--'^SIN(I''-ARES)  +  (MESH( 

C  FN,6)+(FLOAT(STEP)-0.  25)"(MESH(SN,6) -MESH(FN,6) )/ 

C  FLOAT(STEP)-YORG)'''COS(I''^ARES) 

ELSE 

PSI  =  PSISP(FN)+FLOAT(II-l)'nPSISP(SN)-PSISP(FN))/ 
C  FLOAT(STEP) 

DPSI  =  DPSISP(FN)+FL0AT(II-1)'HDPSISP(SN)-DPSISP 
C  (FN))/FLOAT(STEP) 
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DOT  =  N0RM11'>SIN(I'''ARES)  +  N0RM12'''C0S(  I*ARES) 
DOTl  =  (MESH(FN,5)+FLOAT(II-l)''-(MESH(SN,5)-MESH(FN, 
C  5))/FL0AT(STEP)-X0RG)'-SIN(I*ARES)  +  (MESH(FN,6)  + 

C  FLOAT(II-l)^HMESH(SN,6)-MESH(FN,6))/FLOAT(STEP)- 

C  YORG)'>COS(I'^ARES) 

ENDIF 

IF((II.EQ.  1).0R.  (II.  EQ.  STEP+1))  THEN 

TEMP=(J-'^DPSI+D0T*PSI)*(EXP(J'VD0Tl))''^DIST/(2.0* 
C  STEP) 

ELSE 

TEMP=(J*DPSI+D0T'VPSI)*(EXP(J*D0T1))*DIST/ 
C  (STEP) 

ENDIF 

INTEGRAL  =  INTEGRAL  +  TEMP 
15  CONTINUE 

20         CONTINUE 

SIGMA  =  ((CABS(INTEGRAL))**2.0)/(4.  0*K) 
WRITE(50,1010)  I+l,  (I^'-ARES-nSO.O/PI),  SIGMA 
30    CONTINUE 
C 

VRITE(6,''0  '  ' 

WRITE(6,-"-)  '  <«  FIELD  PATTERN  STORED  IN  FFPAT.DAT  >»  ' 
WRITE(6,''^)  '  ' 
C 

1000  FORMAT( IX, 'INTEGRAL  ',13,',  OUT  OF  ',13,'  COMPLETED') 
1010  F0RMAT(1X,I3,2X,F6.  2,2X,E14.  8) 
C 

RETURN 
END 
C 

c 

SUBROUTINE  HAN1(Z,H0,H1) 
C  ^ 

C     Computing  Hankel  Functions  for  n=0,l  with 
C     Complex  Argument,  Z.    Direct  Power  Series  Method  for 
C     CABS(Z)  .LE.  5   and  Hankel' s  Asymptotic  Formula  for 
C     CABS(Z)  .GT.  5.    Written  11/6/87  by  M.  A.  Morgan 
C 

INTEGER  M,M2 

REAL  C(34),DM,F(34),G0,P(34),Pi,P2 

COMPLEX  Z,Z2,Z3,Z4,J0,J1,Y0,Y1,AM,CL,P0,P1,Q0,Q1 

COMPLEX  E0,El,X0,Xl,H0,Hl,j 

PI=3. 1415927 

P2=2.  0/PI 

j=(0.  ,1.) 

IF(CABS(Z).LE.5.0)  THEN 
C 

C  Direct  Power  Series  Method 

C 

G0=  1.  78072 
Z2=0.  5''-Z 
CL=CL0G(G0*Z2) 
C 

C   Computing  F(m)  =  m  !   and  P(m)  =  1  +  1/2  +  1/3  +  +  1/m 

C 

F(l)=1.0 
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P(l)  =  1.0 

DO  11  M=2,34 

F(M)=M''^F(M-1) 
P(M)=P(M-1)+1. 0/M 
11  CONTINUE 

C 

C   Computing  Power  Series  Coefficients 
C 

DM=-1.0 

DO  22  M=l,34 

C(M)=DM/(F(M)*F(M)) 
DM=-DM 
22  CONTINUE 

C 

C   Computing  JO  and  Jl 
C 

J0=(1.  ,0.  ) 

J1=(0.  ,0.  ) 

M=0 
33  M=M+1 

M2=2'^M 

AM=C(M)'V(Z2**M2) 

JO=JO+AM 

J1=J1-M'>AM 

IF((CABS(AM).GT.  l.OE-10).  AND.  (M.  LT.  34))  GO  TO  33 

J1=J1/Z2 
C 

C   Computing  YO  and  Yl 
C 

M=0 

YO=CL*JO 

Y1=Z2"CL"J1-0.5'^J0 
44  M=M+1 

M2=2"M 

AM=C(M)''^P(M)'nZ2'^''^M2) 

YO=YO-AM 

Y1=Y1+M-''AM 

IF((CABS(AM).GT.  1.  OE-10).  AND.  (M.  LT.  34))  GO  TO  44 

Y0=P2''<-Y0 

Y1=P2^VY1/Z2 

HO=JO-j-'-YO 

Hl=Jl-j-'-Yl 

RETURN 
ELSE 
C 

C   Hankel'  Asymptotic  Formula  (Abram.  &  Stegun  p.  364 
C 

Z2=Z^-Z 

Z3=Z'>Z2 

Z4=Z"Z3 

P0=1. 0-. 0703125/Z2+.  1121521/Z4 

Q0=-. 125/Z+.  0732422/Z3 

Pl=l. 0+. 1171875/Z2-.  1441956/Z4 

Ql=. 375/Z-. 10253906/Z3 

XO=(Z-.  25'^PI) 

X1=(Z-. 75-PI) 
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EO=CEXP(-j->XO) 

El=CEXP(-j"Xl) 

AM=CSQRT(P2/Z) 

HO=AM-'-(PO-j'--QO)'^EO 

Hl=AM"(Pl-j-'--Ql)->El 
ENDIF 
C 

RETURN 
END 
C 

c 

SUBROUTINE  CSMINV(A,NDIM,N,DETERM,COND, INORM) 
C 

C     INORM  -  FLAG  TO  NORMALIZE  COLUMNS  AND  ROWS  OF  MATRIX  A 
C     MATRIX  NORMALIZATION  BY  M.  A.  MORGAN 
C     APRIL  24,1978 
C 

C     A      -  MATRIX  TO  INVERT  -  INPUT/ OUTPUT 

C     NDIM   -  -  INPUT 

C     N      -  -  INPUT 

C     DETERM  -  DETERMINATE  OF  A  -  OUTPUT 

C     COND    -  CONDITION  NUMBER  OF  A       -  OUTPUT 
C     INORM   -  INTEGER  NORMALIZATION  FLAG   -  INPUT 
C 
C 

COMPLEX  A(100, 100), PIVOT(100),AMAX,T, SWAP, DETERM, U 
INTEGER  I,J,K,L,IPIVOT(100),INDEX(100,2),IROW,ICOLUM,L1,JROW 
INTEGER  JCOLUM,N, INORM 

REAL  TEMP, ALPHA( 100 ),COL( 100 ),ROW( 100 ),AJK,SUMAXA,SUMROW,SUMAXI 
C 

IF(NDIM. GT. 100)  THEN 

WRITEC^S^-)  '  ERROR  IN  INVERTION  CALL...  DIMENSION  >  100  ' 
STOP 
ENDIF 
IF(N. GT. NDIM)  THEN 

WRITE(--'s^-)  •  ERROR  IN  INVERTION  CALL...  N  >  MAX  DIM.  ' 
STOP 
ENDIF 

IF( INORM.  NE.  1)  GO  TO  7 
DO-  3  K  =  1 ,  N 

COL(K)  =  0.  0 
DO  1  J  =  1,N 

AJK  =  CABS(A(J,K)) 

IF(AJK.GT.  COL(K))  COL(K)  =  AJK 

1  CONTINUE 

DO  2  J  =  1,  N 

2  A(J,K)  =  A(J,K)/COL(K) 

3  CONTINUE 

C     ROW  NORMALIZING 
DO  6  J  =  1,  N 

ROW(J)  =0.0 
DO  4  K  =  1,  N 

AJK  =  CABS(A(J,K)) 

IF(AJK. GT. ROW(J))  ROW(J)  =  AJK 

4  CONTINUE 
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DO  5  K  =  1,  N 

5  A(J,K)  =  A(J,K)/ROW(J) 

6  CONTINUE 

7  CONTINUE 

DETERM  =  CMPLX(1.  0,0.  0) 
SUMAXA  =  0. 0 
DO  20  J  =  1,  N 

ALPHA(J)  =0.0 
SUMROW  =  0. 0 
DO  10  I  =  1,  N 

ALPHA(J)  =  ALPHA(J)  +  ACJ,!)''^  CONJG(  A(  J,  I) ) 
10         SUMROW  =  SUMROW  +  CABS (A( J, I)) 
ALPHA(J)  =  SQRT(ALPHA(J)) 
IF(SUMROW.GT. SUMAXA)  SUMAXA  =  SUMROW 
20    IPIVOT(J)  =  0 

DO  600  I  =  1,  N 
C 
C 

AMAX  =  CMPLX(0. 0,0.  0) 
DO  105  J  =  1,  N 

IF  (IPIVOT(J)-l)  60,  105,  60 
60  DO  100  K  =  1,  N 

IF  (IPIVOT(K)-l)  80,  100,  740 
80  TEMP  =  AMAX''-CONJG(AMAX)  -  A(  J,K)*CONJG(  A(  J,K) ) 

IF(TEMP)85,85,100 
85  IROW  =  J 

ICOLUM  =  K 
AMAX  =  A(J,K) 
100  CONTINUE 

105        CONTINUE 

IPIVOTC ICOLUM)  =  IPIVOTC ICOLUM)  +  1 
C 

c 

IF  (IROW-ICOLUM)  140,  260,  140 
140  DETERM  =  -DETERM 

DO  200  L  =  1,  N 

SWAP  =  A(IROW,L) 
A(IROW,L)  =  A( ICOLUM, L) 
200  A( ICOLUM, L)  =  SWAP 

SWAP  =  ALPHA(IROW) 

ALPHA(IROW)  =  ALPHA( ICOLUM) 

ALPHA( ICOLUM)  =  SWAP 
260  INDEX(I,1)  =  IROW 

INDEX(I,2)  =  ICOLUM 

PIVOT(I)  =  A( ICOLUM, ICOLUM) 

U  =  PIVOT(I) 

DETERM  =  DETERM''^U 

DETERM  =  DETERM/ALPHA( ICOLUM) 

TEMP  =  PIVOT(I)''^CONJG(PIVOT(I)) 

IF(TEMP)330,720,330 
C 
C 
330  A( ICOLUM, ICOLUM)  =  CMPLX( 1. 0,0. 0) 

DO  350  L  =  1,  N 

U  =  PIVOT(I) 
350  A( ICOLUM, L)  =  A( ICOLUM, L)/U 
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c 

c 

380  DO  550  LI  =  1,  N 

IF(Ll-ICOLUM)  400,  550,  400 
400  T  =  A(L1,IC0LUM) 

A(L1,IC0LUM)  =  CMPLX(  0.0,0.0) 

DO  450  L  =  1,  N 

U  =  A(ICOLUM,L) 
450  A(L1,L)  =  A(L1,L)  -  V*T 

550        CONTINUE 
600   CONTINUE 
C 

c 

620   DO  710  I  =  1,  N 

L  =  N  +  1  -  I 

IF  (INDEX(L,1)  -  INDEX(L,2))  630,  710,  630 
630  JROW  =  INDEX(L,1) 

JCOLUM  =  INDEX(L,2) 
DO  705  K  =  1,  N 

SWAP  =  A(K,JROW) 
A(K,JROW)  =  A(K, JCOLUM) 
A(K, JCOLUM)  =  SWAP 
705  CONTINUE 

710    CONTINUE 

SUMAXI  =  0. 0 
DO  910  I  =  1,  N 

SUMROW  =  0.  0 
DO  900  J  =  1,  N 
900        SUMROW  =  SUMROW  +  CABS(A(I,J)) 

IF(SUMROW.GT. SUMAXI)  SUMAXI  =  SUMROW 
910   CONTINUE 

COND  =  SUMAXA'VSUMAXI 
IFdNORM.  NE.  1)  GO  T6~^5^ 
DO  950  K  =  1,  N 

DO  950  J  =  1,  N 
950        A(J,K)  =  A(J,K)/(ROW(K)''^COL(J)) 
955    CONTINUE 

RETURN 
720   WRITE(^'S730) 

730   FORMATC  '  MATRIX  IS  SINGULAR') 
740   RETURN 

END 
C 

c 
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APPENDIX  F.     DIELECTRIC  CYLINDER  SCATTERING  PROGRAM 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C     PLANEWAVE  SCATTERING  BY  A  DIELECTRIC  CYLINDER  C 

C     E  -  WAVE   (TM  CASE)  C 

C     H  -  WAVE   (TE  CASE)  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

COMPLEX  GAMMA(0: 200,2) ,SIGMA(2) ,ER,MU,KR,YR,ZR 

C0MPLEX''-16  JA(0:  200),DJA(0:  200),KRRA 

COMPLEX^ne  J(0:  200),DJ(0:  200) 

REAL---8  Y(0:  200)  ,DY(0:  200),YA(0:  200),DYA(0:  200)  ,JB(0:  200) 

REAL^'-8  DJB(0:  200)  ,RA,KORA 

REAL  K0,PI,SIGMAN(2),A,B 

INTEGER  I, I I, MODE, ARES, N 
C 

0PEN(3,FILE='C:  MSFORT  DECTEM.DAT') 

PI  =  3. 1415927 
C 

WRITEC-'-,^-)  'PERMITTIVITY  FORMAT  IS   "  a  +  jb,  "  ' 

WRITE( •-,•-)  'Enter  Dielectric  Constant  (REAL  PART,  a)  ' 

READ(^'-,^>)  A 

WRITEC^S-'-)  'Enter  Dielectric  Constant  (IMAGINARY  PART,  b)  ' 

READ(-'S^'-)  B 

ER  =  CMPLX(A,B) 
C 

WRITE (^-,^'>-)  'PERMEABILITY  FORMAT  IS   "  a  +  jb,  "  ' 

WRITEO'',")  'Enter  Permeability  Constant  (REAL  PART,  a)  ' 

RE  AD  (-•'-,'>)  A 

WRITE(''-,--'0  'Enter  Permeability  Constant  (IMAGINARY  PART,  b)  ' 

READ(-'-,^'0  B 

MU  =  CMPLX(A,B) 
C 

WRITE(''-,''0  'Enter  the  wave  number  (ko)  ' 

READ(^'-,^'-)  KO 
C 

WRITE(^',^)  'Enter  Cylinder  Radius  (IN  WAVENUMBER  UNITS)   ' 

WRITE(''S''0  'WARNING:  Do  Not  Enter  Zero  !  ' 

RE  AD  (-•'-,''-)  KORA 
C 

KR  =  CSQRT(MU''^ER) 

YR  =  CSQRT(ER/MU) 

ZR  =  CSQRT(MU/ER) 

RA  =  KORA/KO 

KRRA  =  KR^^RA 
C 

WRITE(*,'>)  'Enter  No.  of  Modes:    ' 

READ(''s*)  MODE 

WRITE(",")  'Enter  the  angular  resolution 

READ(-'^''0  ARES 
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WRITE(3,*)  'Cylinder  Scattering  vs.  angle' 
WRITE(3,110)  ER,MU,RA,MODE,KO,KRRA 
WRITE(3,''0  '  ' 
C 

CALL  BES(MODE,K0RA,JB,Y,DJB,DY) 
CALL  DCBJNS  (DCMPLX(KORA,0. OD+00) ,MODE, J,DJ) 
CALL  DCBJNS  (  (KRRA'^KO)  ,MODE  ,  JA,DJA) 
WRITE(''S''0  '  RETURNED  FROM  FINAL  BESSEL  CALL' 
C 

C     CALCULATING  GAMMAs 
C 

DO  10,  N  =  0,  MODE 

GAMMA(N,1)  =  (JA(N)'>DJ(N)-YR*DJA(N)'VJ(N))/(JA(N)*CMPLX(DJ(N) 
C      ,-DY(N))  -  YR-'-DJA(N)''^CMPLX(J(N),-Y(N))) 

GAMMA(N,2)  =  ( JA(N)''-DJ(N) -ZR'VDJA(N)*J(N)  )/( JA(N)*CMPLX(DJ(N) 
C     ,-DY(N))  -  ZR'''DJA(N)'>CMPLX(J(N),-Y(N))) 
WRITE (>v,  1000)  N,GAMMA(N,1),GAMMA(N,2) 
10    CONTINUE 
C 

C     CALCULATING  SIGMAs 
C 

DO  30,  I  =  180,  -180,  -ARES 
DO  40,  II  =  1,  2 

SIGMA(II)  =  (0. 0,0. 0) 
DO  20,  N  =  1,  MODE 

SIGMA(II)  =  SIGMA(II)+2.  0'^GAMMA(N,II)*C0S(N'^PI''^I/180.  0) 
20         CONTINUE 

SIGMA(II)  =  SIGMA(II)  +  GAMMA(0,II) 
SIGMAN(II)  =  ((4.  0/K0)-HCABS(SIGMA(II)))'>*2) 
40      CONTINUE 

WRITE(3,120)  I,  SIGMAN(l),  SIGMAN(2) 
30    CONTINUE 
C 

c 

110   FORMAT(lX,'Er  =         ' ,F6. 4 , 1X,F6.  4, / , 

C  Mu  =         '  ,F6.4,1X,F6.4,/, 

C'  RADIUS  (METERS)  =         ',F8.5,/, 

C'  MAX  MODE        =         ' ,13,/, 

C  KO  =         ',F8.5,/, 

C'  KrRa  =         ' ,F8. 5 , 1X,F8. 5) 

120   F0RMAT(1X,I4,2(3X,E14.4)) 
1000   F0RMAT(1X,I3,1X,2(E12.4,1X,E12.4,4X)) 
C 
C 

STOP 
END 
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APPENDIX  G.     SOFTWARE  SOURCES 


1.  DISSPLA 

Integrated  Software  Systems  Corporation 
10505  Sorrento  Vallev  Road 
San  Diego,  CA  92121 

2.  CURVE-DIGITIZER 
West  Coast  Consultants 

4202  Genesee  Avenue.  Suite  309 
San  Diego,  CA  92117 

3.  Microsoft  FORTR.'\N 
16011  NE  36th  Wav 
BOX  97017 
Redmond,  WA  98073 

4.  Microwav  NDP  FORTRAN 
FOB  79 

Kingston.  MA  02364 

5.  Prof  M.A.  Morgan.  Code  62Mw 
Naval  Postgraduate  School 
Monterey,  CA  93943 

6.  LTT.B.  Welch  III 
c  0  T.B.  Welch  Jr. 
1318  Walthour  Road 
Savannah,  GA  31410 
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